#include "Riostream.h"
#include <cstdlib>
#include <iostream>
#include <iomanip>
#include <string>

#include "TFile.h"
#include "TTree.h"
#include "TChain.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "TPad.h"
#include "TStyle.h"
#include "TROOT.h"
#include "TMath.h"
#include "TRandom.h"
#include "TFormula.h"
#include "TSystem.h"
#include "TEnv.h"

#include "CommandLineInterface.hh"
#include "Barrel.hh"
#include "Annular.hh"
#include "Germanium.hh"
#include "SRIMloader.hh"

using namespace TMath;
using namespace std;

fstream Nfile;

ClassImp(Barrel);
ClassImp(Annular);
ClassImp(Germanium);

SRIM srim;

int main(int argc, char* argv[]){
	vector<char*> InputFiles;
	vector<char*> Histos;
	char* OutputFile = NULL;
	char* MBFile = NULL;
	char* cutfile = NULL;
	char* srimdir = NULL;
	char* SetFile = NULL;
	bool verbose = false;
	int tiltcorr = 0;
	int sim = 0;
	
	const Double_t pi = TMath::Pi();

	Int_t Np, Nd, Nt, Np_FB, Nd_FB, Nt_FB, Np_CD, Nd_CD, Nt_CD, tBe;
	Int_t par[3][5], ID[3][5];
	Int_t Apart, Atarget;
	Float_t Ein, thetain, phiin, xin, yin;
	Float_t Pxin, Pyin, Pzin, Pin;
	Float_t ttarget;
	Float_t theta, phi, thetadet, theta_cm[3][5], theta_coin;
	Float_t thetaN, thetap[3][5], phip[3][5], thetain_coin, phiin_coin;
	Float_t x, y, z, r, alpha;
	Float_t xpad[4], ypad[4], zpad[4];
	Float_t Px, Py, Pz, Pp, PN;
	Float_t vcm, vout[2], gamma[3][5], beta[3][5];
	Float_t Eex[3][5], Edet[3], Efor[3][5], EN, Q, Eex_coin, Ein_coin;
	Float_t EGamDS, EGamDS_2, EGamma_total, EGamDS_total;
	Float_t dEcor, Ecor, Rafter;
	Float_t m[110], MBTheta[8][3][6];
	Double_t time[3][5];
	Float_t CDThreshold[4];
	Float_t PT;
	Float_t a, b, c;

	// Masses
	m[0] = 1.00866491574*931.494e3;
	m[1] = 1.00782503207*931.494e3;
	m[2] = 2.01410177785*931.494e3;
	m[3] = 3.01604927767*931.494e3;
	m[4] = 4.00260325415*931.494e3;
	m[9] = 9.012182201*931.494e3;
	m[10] = 10.013533818*931.494e3;
	m[11] = 11.021657749*931.494e3;
	m[12] = 12.026920737*931.494e3; // 12Be
	m[20] = 19.999981315*931.494e3; // 20F
	m[21] = 20.999948951*931.494e3; // 21F
	m[22] = 22.002998815*931.494e3; // 22F
	m[23] = 22.994466904*931.494e3; // 23Ne
	m[107] = 106.905096820*931.494e3;

	CommandLineInterface* interface = new CommandLineInterface();

	interface->Add("-i", "inputfiles", &InputFiles);
	interface->Add("-o", "outputfile", &OutputFile);
	interface->Add("-M", "miniball angle", &MBFile);
	interface->Add("-C", "cut functions", &cutfile);
	interface->Add("-S", "SRIM files directory", &srimdir);
	interface->Add("-Set", "Setting file", &SetFile);
	interface->Add("-hi", "Histograms to be produced", &Histos);
	interface->Add("-t", "1 tiltcorr histos, 0 no tiltcorr", &tiltcorr);
	interface->Add("-v", "verbose", &verbose);
	interface->Add("-sim", "Set if you are analyzing a simulation",&sim);


	interface->CheckFlags(argc, argv);

	if(InputFiles.size() == 0 || OutputFile == NULL){
		cerr<<"You have to provide at least one input file and the output file!"<<endl;
		exit(1);
	}
	if(Histos.size() == 0){
		cerr << "No histograms will be filled!"<<endl;
		exit(1);
	}
	if(srimdir == NULL){
		cerr << "SRIM directory is needed!"<<endl;
		exit(1);
	}

	cout<<"input file(s):"<<endl;
	for(unsigned int i=0; i<InputFiles.size(); i++){
		cout<<InputFiles[i]<<endl;
	}

	cout<<"output file: "<<OutputFile<< endl;

	// Loading settings
	TEnv* Settings = new TEnv(SetFile);
	
	// Beam energy
	Ein = Settings->GetValue("E.beam",2.85e3);
	Ein = Ein*11;
	thetain = Settings->GetValue("theta.beam",0.);
	thetain = thetain/180.*pi;
	phiin = Settings->GetValue("phi.beam",0.);
	phiin = phiin/180.*pi;
	xin = Settings->GetValue("x.beam",0);
	yin = Settings->GetValue("y.beam",0.);

	ttarget = Settings->GetValue("t.target",10.);
	Atarget = Settings->GetValue("A.target",2);

	if (Atarget==12) m[12] = 12.000000000*931.494e3; // 12C

	// Pad settings
	for (Int_t ipad=0;ipad<4;ipad++){
		xpad[ipad] = Settings->GetValue(Form("x.pad.%d",ipad),0);
		ypad[ipad] = Settings->GetValue(Form("y.pad.%d",ipad),0);
		zpad[ipad] = Settings->GetValue(Form("z.pad.%d",ipad),0);
	}

	// A for particles stopped in dE
	Apart = Settings->GetValue("A.particle",2);

	// Generating SRIM graphs for energyloss calculations
	srim.Creategraphs(srimdir);
	//printf("R = %f, E = %f \n",srim.RMylar[4]->Eval(3.0),srim.EMylar[4]->Eval(11.0));
	Rafter = srim.RT[11]->Eval(Ein);
	Ein = srim.ET[11]->Eval(Rafter - ttarget/(2.*cos(thetain)));
	
	Pxin = sqrt(2*m[11]*Ein)*sin(thetain)*cos(phiin);
	Pyin = sqrt(2*m[11]*Ein)*sin(thetain)*sin(phiin);
	Pzin = sqrt(2*m[11]*Ein)*cos(thetain);
	
	


	TEnv *Angels = new TEnv(MBFile);
	for (int i=0; i<8; i++){
		for (int j=0; j<3; j++){
			for (int k=0; k<6; k++){
				MBTheta[i][j][k] = Angels->GetValue(Form("Theta.%d.%d.%d",i,j,k),0.0);
				//printf("theta = %f \n",MBTheta[i][j][k]);
			}
		}
	}

	TChain* tr;
	tr = new TChain("caltr");
	for(unsigned int i=0; i<InputFiles.size(); i++){
		tr->Add(InputFiles[i]);
	}

	if(tr == NULL){
		cout << "could not find tree caltr in file " << endl;
		for(unsigned int i=0; i<InputFiles.size(); i++){
			cout<<InputFiles[i]<<endl;
		}
		return 3;
	}

	vector<Barrel> *FBarrel = new vector<Barrel>;
	vector<Barrel> *BBarrel = new vector<Barrel>;
	vector<Annular> *FCD = new vector<Annular>;
	vector<Germanium> *Miniball = new vector<Germanium>;
	long long EbisTime;
	long long T1Time;
	long long SuperCycleTime;

	tr->SetBranchAddress("ForwardBarrel",&FBarrel);
	tr->SetBranchAddress("BackwardBarrel",&BBarrel);
	tr->SetBranchAddress("ForwardCD",&FCD);
	tr->SetBranchAddress("Miniball",&Miniball);
	tr->SetBranchAddress("EbisTime",&EbisTime);
	tr->SetBranchAddress("T1Time",&T1Time);
	tr->SetBranchAddress("SuperCycleTime",&SuperCycleTime);

	TFile* outfile = new TFile(OutputFile,"recreate");

	if(outfile->IsZombie()){
		return 4;
	}

	CDThreshold[0] = 500.;
	CDThreshold[1] = 1300.;
	CDThreshold[2] = 500.;
	CDThreshold[3] = 800.;

	// Generating cut functions
	TFormula* For_cut[4][16];
	TFormula* Back_cut[4];
	TFormula* CD_cut[6][4];
	TFormula* CD_PT_cut[6][4];
	TFormula* Be_cut[4];

	TEnv *Cutfunc = new TEnv(cutfile);

	// Forward Barrel
	for (Int_t istrip=0;istrip<16;istrip++){
		for (Int_t icut=0;icut<4;icut++) For_cut[icut][istrip] = new TFormula(Form("For_cut_%d_%d",icut,istrip),Cutfunc->GetValue(Form("For.%d.%d",icut,istrip),"0*x"));
	}
	// 11 Beryllium
	Be_cut[0] = new TFormula("Be_cut_0",Cutfunc->GetValue("For.11","0*x"));
	Be_cut[1] = new TFormula("Be_cut_1",Cutfunc->GetValue("For.12","0*x"));

	// Backward Barrel
	for (Int_t icut=0;icut<4;icut++) Back_cut[icut] = new TFormula(Form("Back_cut_%d",icut),Cutfunc->GetValue(Form("Back.%d",icut),"0*x"));

	// CD
	for (Int_t idet=0;idet<4;idet++){
		for (Int_t icut=0;icut<6;icut++){
			CD_cut[icut][idet] = new TFormula(Form("CD_cut_%d_%d",icut,idet),Cutfunc->GetValue(Form("CD.%d.%d",icut,idet),"0*x"));
			CD_PT_cut[icut][idet] = new TFormula(Form("CD_PT_cut_%d_%d",icut,idet),Cutfunc->GetValue(Form("CD_PT.%d.%d",icut,idet),"0*x"));
		}
	}
	// 11 Beryllium
	Be_cut[2] = new TFormula("Be_cut_2",Cutfunc->GetValue("CD.11","0*x"));
	Be_cut[3] = new TFormula("Be_cut_3",Cutfunc->GetValue("CD.12","0*x"));

	// Prepare for vectorial histograms
//	for (Int_t ih=0;ih<Histos.size();ih++){
		// Barrel
	
		// Position
		TH1F* ForwardBarrel_Pos[4][16];
		TH2F* ForwardBarrel_PosvsStrip[4];
		TH1F* BackwardBarrel_Pos[4][16];
		TH2F* BackwardBarrel_PosvsStrip[4];
	
		// dE
		TH1F* ForwardBarrel_dE[4][16];
		TH1F* ForwardBarrel_dEtot[4];
		TH2F* ForwardBarrel_dEvsStrip[4];
		TH2F* ForwardBarrel_dEvsDetangle[4];
		TH1F* ForwardBarrel_dEangle[4][20];
		TH1F* ForwardBarrel_dEforPos[4][16][16];
		TH2F* ForwardBarrel_dEvsPos[4][16];
		TH1F* BackwardBarrel_dE[4][16];
		TH1F* BackwardBarrel_dEtot[4];
		TH2F* BackwardBarrel_dEvsStrip[4];
		TH1F* BackwardBarrel_dEforPos[4][16];
		TH2F* BackwardBarrel_dEvsPos[4];
	
		// Efor
		TH1F* ForwardBarrel_Efor[4];
		TH2F* ForwardBarrel_EforvsDetangle[4];
		TH2F* ForwardBarrel_EforvsPos[4][16];
	
		// E backdetector
		TH1F* ForwardBarrel_Eback[4];
		TH1F* BackwardBarrel_Eback[4];
		TH2F* ForwardBarrel_dEvsEtot[4];
		TH2F* ForwardBarrel_dEvsE[4][16];
		TH2F* ForwardBarrel_dEvsEangle[4][20];
		TH2F* BackwardBarrel_dEvsEtot[4];
		TH2F* BackwardBarrel_dEvsE[4][16];
		TH2F* BackwardBarrel_dEvsEangle[4][20];
		TH2F* ForwardBarrel_dEvsEcortot[4];
		TH2F* ForwardBarrel_dEvsEcor[4][16];
		TH2F* BackwardBarrel_dEvsEcortot[4];
	
	
		// CD
	
		// dE
		TH1F* CD_dEtot[4];
		TH2F* CD_dEvsStrip[4];
		TH2F* CD_dEvsRing[4];
		TH2F* CD_StripvsRing[4];
		TH2F* CD_EStripvsERing[4];
		TH1F* CD_dEforStrip[4][16];
		TH1F* CD_dEforRing[4][16];
	
		// E backdetector
		TH2F* CD_dEvsEtot[4];
		TH2F* CD_dEvsEcortot[4];
		TH1F* CD_Eback[4];
		TH2F* CD_dEvsEStrip[4][16];
		TH2F* CD_dEvsERing[4][16];
	
	
		// Gammas
	
		// Energy
		TH1F* Gamma_ECrystals1[8][3];
		TH1F* Gamma_ECrystals2[8][3];
		TH1F* Gamma_ESegments[8][3][6];
	
		// Particles
		TH2F* Proton_EtotvsLabangle[4];
		TH2F* Deuteron_EtotvsLabangle[4];
		TH2F* Triton_EtotvsLabangle[4];
		TH1F* Proton_Eex_ForwardBarrel[4];
		TH1F* Deuteron_Eex_ForwardBarrel[4];
		TH1F* Triton_Eex_ForwardBarrel[4];
		TH1F* Proton_Eex_CD[4];
		TH1F* Deuteron_Eex_CD[4];
		TH1F* Triton_Eex_CD[4];
		TH1F* NPT_Eex_BackwardBarrel[4];
		TH1F* Proton_Eex_Angle[20];
		TH1F* Proton_Eex_511_Angle[20];
		TH1F* Proton_Eex_511_Angle_Background[20];
		TH1F* Proton_Eex_2100_Angle[20];
		TH1F* Proton_Eex_2100_Angle_Background[20];
		TH1F* Proton_Eex_2700_Angle[20];
		TH1F* Deuteron_Eex_Angle[20];
		TH1F* Deuteron_Eex_320_Angle[20];
		TH1F* Deuteron_Eex_320_Angle_Background[20];
		TH1F* Triton_Eex_Angle[20];
		TH1F* Triton_Eex_3300_Angle[20];
		TH1F* Triton_Eex_3300_comp_Angle[20];
		TH1F* Triton_Eex_2900_Angle[20];
		TH1F* Triton_Eex_2600_Angle[20];
		TH1F* Triton_Eex_6000_Angle[20];
		TH1F* Triton_Eex_6000_comp_Angle[20];
		TH1F* NPT_Eex_Angle[20];
		TH1F* NPT_Eex_Angle_Back[20];
		TH1F* NPT_Eex_320_Angle[20];
		TH1F* NPT_Eex_320_Angle_Background[20];
		TH1F* NPT_Eex_2100_Angle[20];
		TH1F* NPT_Eex_2100_Angle_Background[20];
		TH1F* NPT_Eex_2600_Angle[20];
		TH1F* NPT_Eex_2700_Angle[20];
		TH1F* NPT_Eex_2900_Angle[20];
		TH1F* NPT_Eex_3300_Angle[20];
		TH1F* NPT_Eex_3300_comp_Angle[20];
		TH1F* NPT_Eex_3300_Angle_Background[20];
		TH1F* NPT_Eex_6000_comp_Angle[20];
		TH1F* NPT_Eex_6000_Angle_Background[20];
		TH1F* Beryllium_Eex_Angle_Stop_FB[20];
		TH1F* Beryllium_Eex_Angle_Stop_CD[20];
		TH1F* Proton_Eex_coin[4];
		TH1F* Deuteron_Eex_coin[4];
		TH1F* Triton_Eex_coin[4];
		TH2F* Proton_EexvsLabangle_ForwardBarrel[4];
		TH2F* Deuteron_EexvsLabangle_ForwardBarrel[4];
		TH2F* Triton_EexvsLabangle_ForwardBarrel[4];
		TH2F* Proton_EexvsLabangle_CD[4];
		TH2F* Deuteron_EexvsLabangle_CD[4];
		TH2F* Triton_EexvsLabangle_CD[4];
		TH2F* Proton_EexvsLabangle_coin[4];
		TH2F* Deuteron_EexvsLabangle_coin[4];
		TH2F* Triton_EexvsLabangle_coin[4];
		TH1F* Proton_ThetaIn_coin[4];
		TH1F* Deuteron_ThetaIn_coin[4];
		TH1F* Triton_ThetaIn_coin[4];
		TH1F* Triton_PhiIn_coin[4];
		TH2F* Proton_ThetaInvsEex_coin[4];
		TH2F* Deuteron_ThetaInvsEex_coin[4];
		TH2F* Triton_ThetaInvsEex_coin[4];
		TH2F* Proton_PhivsPhi_coin[4];
		TH2F* Deuteron_PhivsPhi_coin[4];
		TH2F* Triton_PhivsPhi_coin[4];
		TH2F* Hitpattern_coin[4];
		TH2F* XY_coin[4];
		TH2F* XY_Beam_coin[4];
	
		TH1F* Proton_ESegments[8][3][6];
		TH1F* Deuteron_ESegments[8][3][6];
	
		// No punchthrough
		TH2F* ForwardBarrel_dEvsLabangle[4];
		TH2F* BackwardBarrel_dEvsLabangle[4];
		TH2F* CD_dEvsLabangle[4];
	
	
		// Initiate histograms
	
		// Single histograms
	
		// Angle histograms
		TH1F* ForwardBarrel_Theta = new TH1F("ForwardBarrel_Theta","Theta in forward barrel",180,0.,90.);
		TH1F* ForwardBarrel_Phi = new TH1F("ForwardBarrel_Phi","Phi in forward barrel",360,0.,360.);
		TH1F* ForwardBarrel_Theta_det = new TH1F("ForwardBarrel_Theta_det","detectorangles for forward barrel",180,0.,90.);
		TH1F* BackwardBarrel_Theta = new TH1F("BackwardBarrel_Theta","Theta in backward barrel",100,0.,90.);
		TH1F* BackwardBarrel_Phi = new TH1F("BackwardBarrel_Phi","Phi in backward barrel",360,0.,360.);
		TH1F* CD_Theta = new TH1F("CD_Theta","Theta in CD",180,0.,90.);
		TH1F* CD_Phi = new TH1F("CD_Phi","Phi in CD",360,0.,360.);
	
		// dE vs E histograms
		TH2F* ForwardBarrel_dEvsEtotal = new TH2F("ForwardBarrel_dEvsEtotal","dE vs Eback for all fordet",500,-1.,24.,500,0.,25.);
		TH2F* ForwardBarrel_dEvsEcortotal = new TH2F("ForwardBarrel_dEvsEcortotal","dE vs Eback corrected for all fordet",600,0.,30.,300,0.,15.);
		TH2F* ForwardBarrel_dEvsEcorp = new TH2F("ForwardBarrel_dEvsEcorp","dE vs Eback corrected protons for all fordet",600,0.,30.,300,0.,15.);
		TH2F* ForwardBarrel_dEvsEcord = new TH2F("ForwardBarrel_dEvsEcord","dE vs Eback corrected deuterons for all fordet",600,0.,30.,300,0.,15.);
		TH2F* ForwardBarrel_dEvsEcort = new TH2F("ForwardBarrel_dEvsEcort","dE vs Eback corrected tritons for all fordet",600,0.,30.,300,0.,15.);
		TH2F* BackwardBarrel_dEvsEtotal = new TH2F("BackwardBarrel_dEvsEtotal","dE vs Eback for all backdet",300,0.,6.,500,0.,10.);
		TH2F* BackwardBarrel_dEvsEcortotal = new TH2F("BackwardBarrel_dEvsEcortotal","dE vs Eback corrected for all backdet",400,0.,8.,500,0.,10.);
		TH2F* BackwardBarrel_dEvsEcorp = new TH2F("BackwardBarrel_dEvsEcorp","dE vs Eback corrected protons for all backdet",400,0.,8.,500,0.,10.);
		TH2F* CD_dEvsEtotal = new TH2F("CD_dEvsEtotal","dE vs E for all CD's",600,0.,15.,1000,0.,25.);
		TH2F* CD_dEvsEcortotal = new TH2F("CD_dEvsEcortotal","Corrected dE vs E for all CD's",600,0.,15.,1000,0.,25.);
		TH2F* CD_dEvsEcorp = new TH2F("CD_dEvsEcorp","dE vs Eback corrected protons for all CD",600,0.,15.,1000,0.,25.);
		TH2F* CD_dEvsEcord = new TH2F("CD_dEvsEcord","dE vs Eback corrected deuterons for all CD",600,0.,15.,1000,0.,25.);
		TH2F* CD_dEvsEcort = new TH2F("CD_dEvsEcort","dE vs Eback corrected tritons for all CD",600,0.,15.,1000,0.,25.);
	
		// Particle histograms
		TH2F* Proton_EtotvsLabangletotal = new TH2F("Proton_EtotvsLabangletotal","E vs labangle for protons",180,0.,90.,500,0.,25.);
		TH2F* Proton_EtotvsCmangletotal = new TH2F("Proton_EtotvsCmangletotal","E vs Cmangle for protons",180,0.,180.,500,0.,25.);
		TH2F* Proton_Back_EtotvsLabangletotal = new TH2F("Proton_Back_EtotvsLabangletotal","E vs labangle for protons in backdet",180,90.,180.,300,0.,15.);
		TH1F* Proton_Eextotal = new TH1F("Proton_Eextotal","Excitation energy for 12Be",1000,-2.5,7.5);
		TH1F* Proton_Eex_cointotal = new TH1F("Proton_Eex_cointotal","Excitation energy for 12Be from coincidence",1000,-2.5,7.5);
		TH1F* Proton_EGamma = new TH1F("Proton_EGamma","Gamma energies in coincidence with protons",800,0,8.);
		TH1F* Proton_EGamma_Total = new TH1F("Proton_EGamma_Total","Gamma energies in coincidence with protons",800,0,8.);
		TH1F* Proton_EGamma_1 = new TH1F("Proton_EGamma_1","Gamma energies in coincidence with protons",800,0,8.);
		TH1F* Proton_EGamma_2 = new TH1F("Proton_EGamma_2","Gamma energies in coincidence with protons",800,0,8.);
		TH1F* Proton_EGamma_3 = new TH1F("Proton_EGamma_3","Gamma energies in coincidence with protons",800,0,8.);
		TH1F* Proton_EDS = new TH1F("Proton_EDS","Doppler shifted gamma energies in coincidence with protons",800,0,8.);
		TH1F* Proton_EDS_Total = new TH1F("Proton_EDS_Total","Doppler shifted gamma energies in coincidence with protons with add-back",800,0,8.);
		TH1F* Proton_EDS_1 = new TH1F("Proton_EDS_1","Doppler shifted gamma energies in coincidence with protons",800,0,8.);
		TH1F* Proton_EDS_2 = new TH1F("Proton_EDS_2","Doppler shifted gamma energies in coincidence with protons",800,0,8.);
		TH1F* Proton_EDS_3 = new TH1F("Proton_EDS_3","Doppler shifted gamma energies in coincidence with protons",800,0,8.);
		TH2F* Proton_EexvsEGamma = new TH2F("Proton_EexvsEGamma","Eex vs Egamma for 12Be",500,-2.5,7.5,800,0.,8.);
		TH2F* Proton_EexvsEDS = new TH2F("Proton_EexvsEDS","Eex vs Egamma for 12Be",500,-2.5,7.5,800,0.,8.);
		TH2F* Proton_EexvsEGamma_1 = new TH2F("Proton_EexvsEGamma_1","Eex vs Egamma for 12Be",500,-2.5,7.5,800,0.,8.);
		TH2F* Proton_EexvsEDS_1 = new TH2F("Proton_EexvsEDS_1","Eex vs Egamma for 12Be",500,-2.5,7.5,800,0.,8.);
		TH2F* Proton_EGammavsEGamma = new TH2F("Proton_EGammavsEGamma","Eex vs Egamma for 10Be",800,0.,8.,800,0,8.);
		TH2F* Proton_EDSvsEDS = new TH2F("Proton_EDSvsEDS","Eex vs Egamma for 10Be",800,0.,8.,800,0,8.);
		TH2F* Deuteron_EtotvsLabangletotal = new TH2F("Deuteron_EtotvsLabangletotal","E vs labangle for deuterons",180,0,90,500,0,25.);
		TH2F* Deuteron_EtotvsCmangletotal = new TH2F("Deuteron_EtotvsCmangletotal","E vs Cmangle for Deuterons",180,0,180.,500,0.,25.);
		TH1F* Deuteron_Eextotal = new TH1F("Deuteron_Eextotal","Excitation energy for 11Be",1000,-2.5,7.5);
		TH1F* Deuteron_Eex_cointotal = new TH1F("Deuteron_Eex_cointotal","Excitation energy for 11Be from coincidence",1000,-2.5,7.5);
		TH1F* Deuteron_EGamma = new TH1F("Deuteron_EGamma","Gamma energies in coincidence with deuterons",800,0,8.);
		TH1F* Deuteron_EGamma_Total = new TH1F("Deuteron_EGamma_Total","Gamma energies in coincidence with deuterons",800,0,8.);
		TH1F* Deuteron_EGamma_1 = new TH1F("Deuteron_EGamma_1","Gamma energies in coincidence with Deuterons",800,0,8.);
		TH1F* Deuteron_EGamma_2 = new TH1F("Deuteron_EGamma_2","Gamma energies in coincidence with Deuterons",800,0,8.);
		TH1F* Deuteron_EGamma_3 = new TH1F("Deuteron_EGamma_3","Gamma energies in coincidence with Deuterons",800,0,8.);
		TH1F* Deuteron_EDS = new TH1F("Deuteron_EDS","Doppler shifted gamma energies in coincidence with deuterons",800,0,8.);
		TH1F* Deuteron_EDS_Total = new TH1F("Deuteron_EDS_Total","Doppler shifted gamma energies in coincidence with deuterons with add-back",800,0,8.);
		TH1F* Deuteron_EDS_1 = new TH1F("Deuteron_EDS_1","Doppler shifted gamma energies in coincidence with deuterons",800,0,8.);
		TH1F* Deuteron_EDS_2 = new TH1F("Deuteron_EDS_2","Doppler shifted gamma energies in coincidence with deuterons",800,0,8.);
		TH1F* Deuteron_EDS_3 = new TH1F("Deuteron_EDS_3","Doppler shifted gamma energies in coincidence with deuterons",800,0,8.);
		TH2F* Deuteron_EexvsEGamma = new TH2F("Deuteron_EexvsEGamma","Eex vs Egamma for 11Be",1000,-2.5,7.5,800,0,8.);
		TH2F* Deuteron_EexvsEDS = new TH2F("Deuteron_EexvsEDS","Eex vs Egamma for 11Be",1000,-2.5,7.5,800,0,8.);
		TH2F* Deuteron_EexvsEGamma_1 = new TH2F("Deuteron_EexvsEGamma_1","Eex vs Egamma for 11Be",1000,-2.5,7.5,800,0,8.);
		TH2F* Deuteron_EexvsEDS_1 = new TH2F("Deuteron_EexvsEDS_1","Eex vs Egamma for 11Be",1000,-2.5,7.5,800,0,8.);
		TH2F* Deuteron_EGammavsEGamma = new TH2F("Deuteron_EGammavsEGamma","Eex vs Egamma for 10Be",800,0.,8.,800,0,8.);
		TH2F* Deuteron_EDSvsEDS = new TH2F("Deuteron_EDSvsEDS","Eex vs Egamma for 10Be",800,0.,8.,800,0,8.);
		TH2F* Triton_EtotvsLabangletotal = new TH2F("Triton_EtotvsLabangletotal","E vs labangle for tritons",180,0,90,500,0,25.);
		TH2F* Triton_EtotvsCmangletotal = new TH2F("Triton_EtotvsCmangletotal","E vs Cmangle for tritons",180,0.,180.,500,0.,25.);
		TH1F* Triton_Eextotal = new TH1F("Triton_Eextotal","Excitation energy for 10Be",1000,-2.5,7.5);
		TH1F* Triton_Eex_cointotal = new TH1F("Triton_Eex_cointotal","Excitation energy for 10Be from coincidence",1000,-2.5,7.5);
		TH1F* Triton_Ein1_cointotal = new TH1F("Triton_Ein1_cointotal","Excitation energy for 10Be from coincidence",1000,0.,40.);
		TH1F* Triton_Ein2_cointotal = new TH1F("Triton_Ein2_cointotal","Excitation energy for 10Be from coincidence",1000,0.,40.);
		TH1F* Triton_EGamma = new TH1F("Triton_EGamma","Gamma energies in coincidence with tritons",800,0,8.);
		TH1F* Triton_EGamma_Total = new TH1F("Triton_EGamma_Total","Gamma energies in coincidence with tritons",800,0,8.);
		TH1F* Triton_EGamma_1 = new TH1F("Triton_EGamma_1","Gamma energies in coincidence with tritons",800,0,8.);
		TH1F* Triton_EGamma_2 = new TH1F("Triton_EGamma_2","Gamma energies in coincidence with tritons",800,0,8.);
		TH1F* Triton_EGamma_3 = new TH1F("Triton_EGamma_3","Gamma energies in coincidence with tritons",800,0,8.);
		TH1F* Triton_EDS = new TH1F("Triton_EDS","Doppler shifted gamma energies in coincidence with tritons",800,0,8.);
		TH1F* Triton_EDS_Total = new TH1F("Triton_EDS_Total","Doppler shifted gamma energies in coincidence with tritons with add-back",800,0,8.);
		TH1F* Triton_EDS_1 = new TH1F("Triton_EDS_1","Doppler shifted gamma energies in coincidence with tritons",800,0,8.);
		TH1F* Triton_EDS_2 = new TH1F("Triton_EDS_2","Doppler shifted gamma energies in coincidence with tritons",800,0,8.);
		TH1F* Triton_EDS_3 = new TH1F("Triton_EDS_3","Doppler shifted gamma energies in coincidence with tritons",800,0,8.);
		TH2F* Triton_EexvsEGamma = new TH2F("Triton_EexvsEGamma","Eex vs Egamma for 10Be",1000,-2.5,7.5,800,0,8.);
		TH2F* Triton_EexvsEDS = new TH2F("Triton_EexvsEDS","Eex vs Egamma for 10Be",1000,-2.5,7.5,800,0,8.);
		TH2F* Triton_EexvsEGamma_1 = new TH2F("Triton_EexvsEGamma_1","Eex vs Egamma for 10Be",1000,-2.5,7.5,800,0,8.);
		TH2F* Triton_EexvsEDS_1 = new TH2F("Triton_EexvsEDS_1","Eex vs Egamma for 10Be",1000,-2.5,7.5,800,0,8.);
		TH2F* Triton_EGammavsEGamma = new TH2F("Triton_EGammavsEGamma","Eex vs Egamma for 10Be",800,0.,8.,800,0,8.);
		TH2F* Triton_EDSvsEDS = new TH2F("Triton_EDSvsEDS","Eex vs Egamma for 10Be",800,0.,8.,800,0,8.);
		TH1F* Triton_EDS_coin_3300 = new TH1F("Triton_EDS_coin_3300","Gammas in coincidence with a triton and a 3300keV gamma",800,0.,8.);
		TH1F* NPT_Eextotal = new TH1F("NPT_Eextotal","Excitation energy for 12Be",1000,-2.5,7.5);
		TH1F* NPT_Eextotal_Back = new TH1F("NPT_Eextotal_Back","Excitation energy for 12Be",1000,-6.5,3.5);
		TH1F* NPT_Eex_cointotal = new TH1F("NPT_Eex_cointotal","Excitation energy for 10Be from coincidence",1000,-2.5,7.5);
		TH2F* NPT_EexvsEex1_cointotal = new TH2F("NPT_EexvsEex1_cointotal","Excitation energy for 10Be from coincidence",200,-6,30,200,-2.5,10.);
		TH2F* NPT_EexvsEex2_cointotal = new TH2F("NPT_EexvsEex2_cointotal","Excitation energy for 10Be from coincidence",200,-60,5,200,-2.5,10.);
		TH2F* NPT_EtotvsLabangletotal_coin = new TH2F("NPT_EtotvsLabangletotal_coin","Excitation energy for 10Be from coincidence",180,0,90,200,0,30);
		//TH1F* Flour_Eextotal = new TH1F("Flour_Eextotal","Excitation energy for protons in backdet",200,-2,5.);
		
		TH1F* Deuteron_pr_Time = new TH1F("Deuteron_pr_Time","Number of deuterons pr. 10 sek",360,0,3600.);
		
//		TH1F* Proton_Eex_Gamma_1400 = new TH1F("Proton_Eex_Gamma_1400","Excitation energy gated on 1.35-1.45 MeV gammas",90,-5.e2,4.e3);
		TH1F* Proton_Eex_Gamma_511 = new TH1F("Proton_Eex_Gamma_511","Excitation energy gated on 0.49-.52 MeV gammas",200,-.5,3.5);
		TH1F* Proton_Eex_Gamma_511_Background = new TH1F("Proton_Eex_Gamma_511_Background","Excitation energy gated on 0.49-.52 MeV gammas",200,-.5,3.5);
		TH1F* Proton_Eex_Gamma_2100 = new TH1F("Proton_Eex_Gamma_2100","Excitation energy gated on 2.07-2.15 MeV gammas",200,-.5,3.5);
		TH1F* Proton_Eex_Gamma_2100_1 = new TH1F("Proton_Eex_Gamma_2100_1","Excitation energy gated on 2.07-2.15 MeV gammas",200,-.5,3.5);
//		TH1F* Proton_Eex_Gamma_2350 = new TH1F("Proton_Eex_Gamma_2350","Excitation energy gated on 2.3-2.4 MeV gammas",90,-5.e2,4.e3);
		TH1F* Proton_Eex_Gamma_2700 = new TH1F("Proton_Eex_Gamma_2700","Excitation energy gated on 2.68-2.8 MeV gammas",200,-.5,3.5);
		TH1F* Proton_Eex_Gamma_2700_1 = new TH1F("Proton_Eex_Gamma_2700_1","Excitation energy gated on 2.68-2.8 MeV gammas",200,-.5,3.5);
		TH1F* Deuteron_Eex_Gamma_320 = new TH1F("Deuteron_Eex_Gamma_320","Excitation energy gated on .3-.34 MeV gammas",200,-2.,2.);
		TH1F* Deuteron_Eex_Gamma_320_1 = new TH1F("Deuteron_Eex_Gamma_320_1","Excitation energy gated on .3-.34 MeV gammas",200,-2.,2.);
		TH1F* Triton_Eex_Gamma_3300 = new TH1F("Triton_Eex_Gamma_3300","Excitation energy gated on 3.3-3.5 MeV gammas",200,0.,8.);
		TH1F* Triton_Eex_Gamma_3300_comp = new TH1F("Triton_Eex_Gamma_3300_comp","Excitation energy gated on 2.2-3.5 MeV gammas",200,0.,8.);
		TH1F* Triton_Eex_Gamma_3300_1 = new TH1F("Triton_Eex_Gamma_3300_1","Excitation energy gated on 3.3-3.5 MeV gammas",200,0.,8.);
		TH1F* Triton_Eex_Gamma_6000_comp = new TH1F("Triton_Eex_Gamma_6000_comp","Excitation energy gated on 4.6-6.3.2 MeV gammas",200,0.,8.);
		TH1F* Triton_Eex_Gamma_6000_comp_1 = new TH1F("Triton_Eex_Gamma_6000_comp_1","Excitation energy gated on 4.6-6.2 MeV gammas",200,0.,8.);

		TH1F* NPT_Eex_Gamma_320 = new TH1F("NPT_Eex_Gamma_320","Excitation energy gated on 0.31-0.33 MeV gammas",500,-2.,8.);
		TH1F* NPT_Eex_Gamma_320_Background = new TH1F("NPT_Eex_Gamma_320_Background","Excitation energy gated on 0.3-0.31 og 0.33-0.34 MeV gammas",500,-2.,8.);
		TH1F* NPT_Eex_Gamma_2100 = new TH1F("NPT_Eex_Gamma_2100","Excitation energy gated on 2.06-2.16 MeV gammas",500,-2.,8.);
		TH1F* NPT_Eex_Gamma_2100_Background = new TH1F("NPT_Eex_Gamma_2100_Background","Excitation energy gated on 2.16-2.26 MeV gammas",500,-2.,8.);
		TH1F* NPT_Eex_Gamma_2600 = new TH1F("NPT_Eex_Gamma_2600","Excitation energy gated on 2.52-2.68 MeV gammas",500,-2.,8.);
		TH1F* NPT_Eex_Gamma_2700 = new TH1F("NPT_Eex_Gamma_2700","Excitation energy gated on 2.69-2.79 MeV gammas",500,-2.,8.);
		TH1F* NPT_Eex_Gamma_2900 = new TH1F("NPT_Eex_Gamma_2900","Excitation energy gated on 2.8-3.0 MeV gammas",500,-2.,8.);
		TH1F* NPT_Eex_Gamma_3300 = new TH1F("NPT_Eex_Gamma_3300","Excitation energy gated on 3.25-3.5 MeV gammas",500,-2.,8.);
		TH1F* NPT_Eex_Gamma_3300_comp = new TH1F("NPT_Eex_Gamma_3300_comp","Excitation energy gated on 2.5-3.5 MeV gammas",500,-2.,8.);
		TH1F* NPT_Eex_Gamma_3300_Background = new TH1F("NPT_Eex_Gamma_3300_Background","Excitation energy gated on 3.5-3.6 MeV gammas",500,-2.,8.);
		TH1F* NPT_Eex_Gamma_6000_comp = new TH1F("NPT_Eex_Gamma_6000_comp","Excitation energy gated on 4.2-6.2 MeV gammas",500,-2.,8.);
		TH1F* NPT_Eex_Gamma_6000_Background = new TH1F("NPT_Eex_Gamma_6000_Background","Excitation energy gated on 6.2-6.6 MeV gammas",500,-2.,8.);

		TH1F* NPT_Eex_coin_Gamma_2700 = new TH1F("NPT_Eex_coin_Gamma_2700","Excitation energy for 12Be",200,0.,30.);

		TH2F* Proton_EtotvsCmangle_Gamma_2100 = new TH2F("Proton_EtotvsCmangle_Gamma_2100","E vs Cmangle for protons gated on 2.07-2.15MeV",180,0.,180.,500,0.,25.);
		TH2F* Proton_EtotvsCmangle_Gamma_2700 = new TH2F("Proton_EtotvsCmangle_Gamma_2700","E vs Cmangle for protons gated on 2.68-2.7MeV",180,0.,180.,500,0.,25.);
		TH2F* Deuteron_EtotvsCmangle_Gamma_320 = new TH2F("Deuteron_EtotvsCmangle_Gamma_320","E vs Cmangle for deuterons gated on 0.3-0.35MeV",180,0.,180.,500,0.,25.);
		TH2F* Triton_EtotvsCmangle_Gamma_3300 = new TH2F("Triton_EtotvsCmangle_Gamma_3300","E vs Cmangle for tritons gated on 3.3-3.5MeV",180,0.,180.,500,0.,25.);
		TH2F* Triton_EtotvsCmangle_Gamma_6000 = new TH2F("Triton_EtotvsCmangle_Gamma_6000","E vs Cmangle for tritons gated on 4.6-6.2MeV",180,0.,180.,500,0.,25.);
		TH2F* NPT_EtotvsCmangle_Gamma_2100 = new TH2F("NPT_EtotvsCmangle_Gamma_2100","E vs Cmangle for protons gated on 2.07-2.15MeV",180,0.,180.,500,0.,25.);
		TH2F* NPT_EtotvsCmangle_Gamma_2700 = new TH2F("NPT_EtotvsCmangle_Gamma_2700","E vs Cmangle for protons gated on 2.68-2.7MeV",180,0.,180.,500,0.,25.);
		TH2F* NPT_EtotvsCmangle_Gamma_320 = new TH2F("NPT_EtotvsCmangle_Gamma_320","E vs Cmangle for deuterons gated on 0.3-0.35MeV",180,0.,180.,500,0.,25.);
		TH2F* NPT_EtotvsLabangle_Gamma_320 = new TH2F("NPT_EtotvsLabangle_Gamma_320","E vs Cmangle for deuterons gated on 0.3-0.35MeV",180,0.,90.,500,0.,10.);
		TH2F* NPT_EtotvsCmangle_Gamma_3300 = new TH2F("NPT_EtotvsCmangle_Gamma_3300","E vs Cmangle for tritons gated on 3.3-3.5MeV",180,0.,180.,500,0.,25.);
		TH2F* NPT_EtotvsCmangle_Gamma_6000 = new TH2F("NPT_EtotvsCmangle_Gamma_6000","E vs Cmangle for tritons gated on 4.6-6.2MeV",180,0.,180.,500,0.,25.);


		TH1F* Proton_MBTheta_2100 = new TH1F("Proton_MBTheta_2100","Angular distribution for 2100 gammas",90,0,180);
		TH1F* Proton_MBTheta_2700 = new TH1F("Proton_MBTheta_2700","Angular distribution for 2700 gammas",90,0,180);
		TH1F* Deuteron_MBTheta_320 = new TH1F("Deuteron_MBTheta_320","Angular distribution for 320 gammas",90,0,180);
		TH1F* Triton_MBTheta_3300 = new TH1F("Triton_MBTheta_3300","Angular distribution for 3300 gammas",90,0,180);

//		TH2F* Proton_EtotvsLabangle_Gamma_1000 = new TH2F("Proton_EtotvsLabangle_Gamma_1000","E vs theta for protons gated on 1000 MeV gammas",180,0.,90.,1.e2,0.,1.e4);
		
		TH2F* Hitpattern_ForwardBarreltotal = new TH2F("Hitpattern_ForwardBarreltotal","(x,y) for coincidence particles",500,-100,100,500,-100,100);
		TH2F* Hitpattern_CDtotal = new TH2F("Hitpattern_CDtotal","(x,y) for coincidence particles",400,-40,40,400,-40,40);
		TH2F* Hitpattern_cointotal = new TH2F("Hitpattern_cointotal","(x,y) for coincidence particles",400,-40,40,400,-40,40);
		TH2F* XY_cointotal = new TH2F("XY_cointotal","(x,y) for coincidence particles",200,-1,1,200,-1,1);
		TH2F* XY_Beam_cointotal = new TH2F("XY_Beam_cointotal","(x,y) for coincidence particles",500,-50,50,500,-50,50);
		TH2I* DetvsDet_coin = new TH2I("DetvsDet_coin","Detector id vs detector id for coincidence events",6,-1,5,6,-1,5);

		// No punchthrough
		TH2F* ForwardBarrel_dEvsLabangletotal = new TH2F("ForwardBarrel_dEvsLabangletotal","E vs labangle in for dE detector",180,0,90,600,0,30.);
		TH2F* ForwardBarrel_dEvsCmangletotal = new TH2F("ForwardBarrel_dEvsCmangletotal","E vs labangle in for dE detector",180,0,180.,600,0,30.);
		TH2F* BackwardBarrel_dEvsLabangletotal = new TH2F("BackwardBarrel_dEvsLabangletotal","E vs labangle in back dE detector",180,90,180,300,0,15.);
		TH2F* CD_dEvsLabangletotal = new TH2F("CD_dEvsLabangletotal","E vs Labangle in CD dE detector",70,0,35,500,0,50.);
		TH2F* All_dEvsLabangletotal = new TH2F("All_dEvsLabangletotal","E vs Labangle in All dE detector",180,0,pi/2.,500,0,50.e3);
		TH1F* NPT_EGamma = new TH1F("NPT_EGamma","Gamma energies in coincidence with particles stopped in dE detectors",800,0,8.);
		TH1F* NPT_EGamma_Total = new TH1F("NPT_EGamma_Total","Gamma energies in coincidence with particles stopped in dE detectors",800,0,8.);
		TH1F* NPT_EGamma_1 = new TH1F("NPT_EGamma_1","Gamma energies in coincidence with particles stopped in dE detectors",800,0,8.);
		TH1F* NPT_EGamma_2 = new TH1F("NPT_EGamma_2","Gamma energies in coincidence with particles stopped in dE detectors",800,0,8.);
		TH1F* NPT_EGamma_3 = new TH1F("NPT_EGamma_3","Gamma energies in coincidence with particles stopped in dE detectors",800,0,8.);
		TH1F* NPT_EDS = new TH1F("NPT_EDS","Doppler shifted gamma energies in coincidence with particles stopped in dE detectors",800,0,8.);
		TH1F* NPT_EDS_Total = new TH1F("NPT_EDS_Total","Doppler shifted gamma energies in coincidence with particles stopped in dE detectors",800,0,8.);
		TH1F* NPT_EDS_1 = new TH1F("NPT_EDS_1","Doppler shifted gamma energies in coincidence with particles stopped in dE detectors",800,0,8.);
		TH1F* NPT_EDS_2 = new TH1F("NPT_EDS_2","Doppler shifted gamma energies in coincidence with particles stopped in dE detectors",800,0,8.);
		TH1F* NPT_EDS_3 = new TH1F("NPT_EDS_3","Doppler shifted gamma energies in coincidence with particles stopped in dE detectors",800,0,8.);
		TH1F* NPT_EGamma_Low = new TH1F("NPT_EGamma_Low","Gamma energies in coincidence with particles stopped in dE detectors",1000,0,2.);
		TH1F* NPT_EDS_Low = new TH1F("NPT_EDS_Low","Doppler shifted gamma energies in coincidence with particles stopped in dE detectors",1000,0,2.);
		TH1F* NPT_EGamma_High = new TH1F("NPT_EGamma_High","Gamma energies in coincidence with particles stopped in dE detectors",1000,2.,8.);
		TH1F* NPT_EDS_High = new TH1F("NPT_EDS_High","Doppler shifted gamma energies in coincidence with particles stopped in dE detectors",1000,2.,8.);
		TH1F* NPT_EGamma_Back = new TH1F("NPT_EGamma_Back","Gamma energies in coincidence with particles stopped in dE detectors",800,0,8.);
		TH1F* NPT_EDS_Back = new TH1F("NPT_EDS_Back","Doppler shifted gamma energies in coincidence with particles stopped in dE detectors",800,0,8.);
		TH2F* NPT_EexvsEGamma = new TH2F("NPT_EexvsEGamma","Excitation energy vs gamma energy",200,-5.,15.,200,0,8.);
		TH2F* NPT_EexvsEDS = new TH2F("NPT_EexvsEDS","Excitation energy vs doppler corrected gamma energy",200,-5.,15.,200,0,8.);
		TH2F* NPT_EexvsEGamma_1 = new TH2F("NPT_EexvsEGamma_1","Excitation energy vs gamma energy",200,-5.,15.,200,0,8.);
		TH2F* NPT_EexvsEDS_1 = new TH2F("NPT_EexvsEDS_1","Excitation energy vs doppler corrected gamma energy",200,-5.,15.,200,0,8.);
		TH2F* NPT_EexvsEGamma_Back = new TH2F("NPT_EexvsEGamma_Back","Excitation energy vs gamma energy",200,-6.5,3.5,200,0,8.);
		TH2F* NPT_EexvsEDS_Back = new TH2F("NPT_EexvsEDS_Back","Excitation energy vs doppler corrected gamma energy",200,-6.5,3.5,200,0,8.);
		TH2F* NPT_EGammavsEGamma = new TH2F("NPT_EGammavsEGamma","Eex vs Egamma for 10Be",800,0.,8.,800,0,8.);
		TH2F* NPT_EDSvsEDS = new TH2F("NPT_EDSvsEDS","Eex vs Egamma for 10Be",800,0.,8.,800,0,8.);

		// 11Be
		TH1F* Be11_Angle = new TH1F("Be11_Angle","Theta for 11Be",180,0,90);
		TH1F* Be11_pr_Time = new TH1F("Be11_pr_Time","Number of 11Be pr. 10 sek",360,0,3600.);
	
		// Gamma
		TH1F* Gamma_Etotal = new TH1F("Gamma_Etotal","Gamma energies",800,0,8.);
		TH1F* Gamma_EDS = new TH1F("Gamma_EDS","Doppler shifted gamma energies",800,0,8.);
	
		// Time
		//
		// Match this with the `Coincidence.Window` and `Coincidence.Reference`
		// parameters in the setting file used for unpacking the events with
		// `TransferEventBuilder`.
		Float_t timediff_diff_min = -500.;
		Float_t timediff_diff_max = 500.;
		// The detector data gives time signals in quanta of 25 ns. The amount of
		// bins should be the amount of 25 ns intervals in the plot. For best
		// results, make the total above specified interval length a multiple of
		// 25.
		Int_t timediff_xbins = TMath::Nint((timediff_diff_max-timediff_diff_min)/25.);
	
		TH2F* Gamma_EvsdtProton = new TH2F("Gamma_EvsdtProton","Gamma energy vs timediff between gamma and proton",timediff_xbins,timediff_diff_min,timediff_diff_max,1e3,0,4.);
		TH2F* Gamma_EvsdtDeuteron = new TH2F("Gamma_EvsdtDeuteron","Gamma energy vs timediff between gamma and deuteron",timediff_xbins,timediff_diff_min,timediff_diff_max,1e3,0,4.);
		TH2F* Gamma_EvsdtTriton = new TH2F("Gamma_EvsdtTriton","Gamma energy vs timediff between gamma and triton",timediff_xbins,timediff_diff_min,timediff_diff_max,1e3,0,4.);

		TH2F* Gamma_EvsdtNPT = new TH2F("Gamma_EvsdtNPT","Gamma energy vs timediff between gamma and NPT particles",timediff_xbins,timediff_diff_min,timediff_diff_max,1e3,0,4.);
		TH2F* Gamma_EvsdtBack = new TH2F("Gamma_EvsdtBack","Gamma energy vs timediff between gamma and particles in back",timediff_xbins,timediff_diff_min,timediff_diff_max,1e3,0,4.);
		TH2F* Gamma_EvsdtCD = new TH2F("Gamma_EvsdtCD","Gamma energy vs timediff between gamma and particles in CD",timediff_xbins,timediff_diff_min,timediff_diff_max,1e3,0,4.);
		TH2F* Gamma_Evsdttot = new TH2F("Gamma_Evsdttot","Gamma energy vs timediff between gamma and particle",timediff_xbins,timediff_diff_min,timediff_diff_max,1e3,0,4.);
		TH1F* Gamma_dtGe = new TH1F("Gamma_dtGe","Decay time of excited Ge state",timediff_xbins,timediff_diff_min,timediff_diff_max);
		TH1F* Gamma_Ecut = new TH1F("Gamma_Ecut","E gamma cut for delayed decay",1e3,2.,2.5);
		
		TH1F* Proton_dt_511 = new TH1F("Proton_dt_511","timediff between gamma and proton for 0p2",timediff_xbins,timediff_diff_min,timediff_diff_max);
	
		// Pad dependent histograms
		for (Int_t i=0;i<4;i++){
			// Position
			ForwardBarrel_PosvsStrip[i] = new TH2F(Form("ForwardBarrel_PosvsStrip_%d",i),Form("Position vs stripnumber for fordet %d",i),16,0,16,220,0,1.1);
			BackwardBarrel_PosvsStrip[i] = new TH2F(Form("BackwardBarrel_PosvsStrip_%d",i),Form("Position vs stripnumber for backdet %d",i),16,0,16,220,0,1.1);
	
			// dE
			ForwardBarrel_dEtot[i] = new TH1F(Form("ForwardBarrel_dEtot_%d",i),Form("Total dE for fordet %d",i),600,0,30.);
			ForwardBarrel_dEvsStrip[i] = new TH2F(Form("ForwardBarrel_dEvsStrip_%d",i),Form("dE vs stripnumber for fordet %d",i),16,0,16,300,0,30.);
			ForwardBarrel_dEvsDetangle[i] = new TH2F(Form("ForwardBarrel_dEvsDetangle_%d",i),Form("dE vs detector angle for fordet %d",i),180,0,pi/2.,300,0,30.);
			BackwardBarrel_dEtot[i] = new TH1F(Form("BackwardBarrel_dEtot_%d",i),Form("Total dE for backdet %d",i),300,0,15.);
			BackwardBarrel_dEvsStrip[i] = new TH2F(Form("BackwardBarrel_dEvsStrip_%d",i),Form("dE vs stripnumber for backdet %d",i),16,0,16,150,0,15.);
			BackwardBarrel_dEvsPos[i] = new TH2F(Form("BackwardBarrel_dEvsPos_%d",i),Form("dE vs position for backdet %d",i),330,0,1.1,250,0,2.5e3);
	
			// Efor
			ForwardBarrel_Efor[i] = new TH1F(Form("ForwardBarrel_Efor_%d",i),Form("Total Efor for fordet %d",i),700,0.,7.);
			ForwardBarrel_EforvsDetangle[i] = new TH2F(Form("ForwardBarrel_EforvsDetangle_%d",i),Form("Efor vs detector angle for fordet %d",i),180,0,pi/2.,500,0.,10.);
	
			// E backdetector
			ForwardBarrel_Eback[i] = new TH1F(Form("ForwardBarrel_Eback_%d",i),Form("Eback for fordet %d",i),600,0.,30.);
			BackwardBarrel_Eback[i] = new TH1F(Form("BackwardBarrel_Eback_%d",i),Form("Eback for Backdet %d",i),300,0.,15.);
			ForwardBarrel_dEvsEtot[i] = new TH2F(Form("ForwardBarrel_dEvsEtot_%d",i),Form("dE vs Eback for fordet %d",i),500,-1.,24.,500,0.,25.);
			BackwardBarrel_dEvsEtot[i] = new TH2F(Form("BackwardBarrel_dEvsEtot_%d",i),Form("dE vs Eback for backdet %d",i),300,0.,6.,500,0.,10.);
			ForwardBarrel_dEvsEcortot[i] = new TH2F(Form("ForwardBarrel_dEvsEcortot_%d",i),Form("dE vs Eback corrected for fordet %d",i),600,0.,30.,300,0.,15.);
			BackwardBarrel_dEvsEcortot[i] = new TH2F(Form("BackwardBarrel_dEvsEcortot_%d",i),Form("dE vs Eback corrected for backdet %d",i),400,0.,8.,500,0.,10.);
	
			// No punchthrough
			ForwardBarrel_dEvsLabangle[i] = new TH2F(Form("ForwardBarrel_dEvsLabangle_%d",i),Form("dE vs labangle for fordet %d",i),180,0,90,600,0,30.);
			BackwardBarrel_dEvsLabangle[i] = new TH2F(Form("BackwardBarrel_dEvsLabangle_%d",i),Form("dE vs labangle for backdet %d",i),180,90,180,300,0,15.);
	
			// particles
			Proton_EtotvsLabangle[i] = new TH2F(Form("Proton_EtotvsLabangle_%d",i),Form("E vs labangle for protons in det %d",i),180,0,90,500,0,25.);
			Deuteron_EtotvsLabangle[i] = new TH2F(Form("Deuteron_EtotvsLabangle_%d",i),Form("E vs labangle for deuterons in det %d",i),180,0,90,500,0,25.);
			Triton_EtotvsLabangle[i] = new TH2F(Form("Triton_EtotvsLabangle_%d",i),Form("E vs labangle for tritons in det %d",i),180,0,90,500,0,25.);
			Proton_Eex_ForwardBarrel[i] = new TH1F(Form("Proton_Eex_ForwardBarrel_%d",i),Form("Excitation energy for protons in fordet %d",i),700,-2,5.);
			Deuteron_Eex_ForwardBarrel[i] = new TH1F(Form("Deuteron_Eex_ForwardBarrel_%d",i),Form("Excitation energy for deuterons in fordet %d",i),500,-2.,3.);
			Triton_Eex_ForwardBarrel[i] = new TH1F(Form("Triton_Eex_ForwardBarrel_%d",i),Form("Excitation energy for tritons in fordet %d",i),1000,-2.5,7.5);
			Proton_EexvsLabangle_ForwardBarrel[i] = new TH2F(Form("Proton_EexvsLabangle_ForwardBarrel_%d",i),Form("Excitation energy for protons in fordet %d",i),180,0,90,700,-2,5.);
			Deuteron_EexvsLabangle_ForwardBarrel[i] = new TH2F(Form("Deuteron_EexvsLabangle_ForwardBarrel_%d",i),Form("Excitation energy for deuterons in fordet %d",i),180,0,90,500,-2.,3.);
			Triton_EexvsLabangle_ForwardBarrel[i] = new TH2F(Form("Triton_EexvsLabangle_ForwardBarrel_%d",i),Form("Excitation energy for tritons in fordet %d",i),180,0,90,1000,-2.5,7.5);
			NPT_Eex_BackwardBarrel[i] = new TH1F(Form("NPT_Eex_BackwardBarrel_%d",i),Form("Excitation energy for particles in backdet %d",i),1000,-8.,2.);
			Proton_Eex_CD[i] = new TH1F(Form("Proton_Eex_CD_%d",i),Form("Excitation energy for protons in fordet %d",i),700,-2,5.);
			Deuteron_Eex_CD[i] = new TH1F(Form("Deuteron_Eex_CD_%d",i),Form("Excitation energy for deuterons in fordet %d",i),500,-2.,3.);
			Triton_Eex_CD[i] = new TH1F(Form("Triton_Eex_CD_%d",i),Form("Excitation energy for tritons in fordet %d",i),1000,-2.5,7.5);
			Proton_EexvsLabangle_CD[i] = new TH2F(Form("Proton_EexvsLabangle_CD_%d",i),Form("Excitation energy for protons in fordet %d",i),180,0,90,700,-2,5.);
			Deuteron_EexvsLabangle_CD[i] = new TH2F(Form("Deuteron_EexvsLabangle_CD_%d",i),Form("Excitation energy for deuterons in fordet %d",i),180,0,90,500,-2.,3.);
			Triton_EexvsLabangle_CD[i] = new TH2F(Form("Triton_EexvsLabangle_CD_%d",i),Form("Excitation energy for tritons in fordet %d",i),180,0,90,1000,-2.5,7.5);
			Proton_Eex_coin[i] = new TH1F(Form("Proton_Eex_coin_%d",i),Form("Excitation energy for protons in fordet %d from coincidence",i),700,-2,5.);
			Deuteron_Eex_coin[i] = new TH1F(Form("Deuteron_Eex_coin_%d",i),Form("Excitation energy for deuterons in fordet %d from coincidence",i),500,-2.,3.);
			Triton_Eex_coin[i] = new TH1F(Form("Triton_Eex_coin_%d",i),Form("Excitation energy for tritons in fordet %d from coincidence",i),1000,-2.5,7.5);
			Proton_EexvsLabangle_coin[i] = new TH2F(Form("Proton_EexvsLabangle_coin_%d",i),Form("Excitation energy for protons in fordet %d",i),180,0,90,700,-2,5.);
			Deuteron_EexvsLabangle_coin[i] = new TH2F(Form("Deuteron_EexvsLabangle_coin_%d",i),Form("Excitation energy for deuterons in fordet %d",i),180,0,90,500,-2.,3.);
			Triton_EexvsLabangle_coin[i] = new TH2F(Form("Triton_EexvsLabangle_coin_%d",i),Form("Excitation energy for tritons in fordet %d",i),180,0,90,1000,-2.5,7.5);
			Proton_ThetaIn_coin[i] = new TH1F(Form("Proton_ThetaIn_coin_%d",i),Form("Incoming angle for for det %d calculated from coincidence",i),160,-20,20);
			Deuteron_ThetaIn_coin[i] = new TH1F(Form("Deuteron_ThetaIn_coin_%d",i),Form("Incoming angle for for det %d calculated from coincidence",i),160,-20,20);
			Triton_ThetaIn_coin[i] = new TH1F(Form("Triton_ThetaIn_coin_%d",i),Form("Incoming angle for for det %d calculated from coincidence",i),160,-20,20);
			Triton_PhiIn_coin[i] = new TH1F(Form("Triton_PhiIn_coin_%d",i),Form("Incoming angle for for det %d calculated from coincidence",i),360,0,360);
			Proton_ThetaInvsEex_coin[i] = new TH2F(Form("Proton_ThetaInvsEex_coin_%d",i),Form("Incoming angle for for det %d calculated from coincidence",i),160,-20,20,700,-2,5.);
			Deuteron_ThetaInvsEex_coin[i] = new TH2F(Form("Deuteron_ThetaInvsEex_coin_%d",i),Form("Incoming angle for for det %d calculated from coincidence",i),160,-20,20,500,-2,3.);
			Triton_ThetaInvsEex_coin[i] = new TH2F(Form("Triton_ThetaInvsEex_coin_%d",i),Form("Incoming angle for for det %d calculated from coincidence",i),160,-20,20,1000,-2.5,7.5);
			Proton_PhivsPhi_coin[i] = new TH2F(Form("Proton_PhivsPhi_%d",i),Form("phi vs phi for det %d",i),360,0,360,360,0,360);
			Deuteron_PhivsPhi_coin[i] = new TH2F(Form("Deuteron_PhivsPhi_%d",i),Form("phi vs phi for det %d",i),360,0,360,360,0,360);
			Triton_PhivsPhi_coin[i] = new TH2F(Form("Triton_PhivsPhi_%d",i),Form("phi vs phi for det %d",i),360,0,360,360,0,360);
	
			Hitpattern_coin[i] = new TH2F(Form("Hitpattern_coin_%d",i),Form("(x,y) for coincidence particles in for det %d",i),400,-40,40,400,-40,40);
			XY_coin[i] = new TH2F(Form("XY_coin_%d",i),Form("(x,y) for coincidence particles in for det %d",i),200,-1,1,200,-1,1);
			XY_Beam_coin[i] = new TH2F(Form("XY_Beam_coin_%d",i),Form("(x,y) for coincidence particles in for det %d",i),200,-20,20,200,-20,20);
	
			// Strip dependent histograms
			for (Int_t j=0;j<16;j++){
				// Position
				ForwardBarrel_Pos[i][j] = new TH1F(Form("ForwardBarrel_Pos_%d_%d",i,j),Form("Position for strip %d in fordet %d",j,i),750.,0.,2.);
				BackwardBarrel_Pos[i][j] = new TH1F(Form("BackwardBarrel_Pos_%d_%d",i,j),Form("Position for strip %d in backdet %d",j,i),750,0.,1.5);
	
				// dE
				ForwardBarrel_dE[i][j] = new TH1F(Form("ForwardBarrel_dE_%d_%d",i,j),Form("dE for strip %d in fordet %d",j,i),600,0,30.);
				ForwardBarrel_dEvsPos[i][j] = new TH2F(Form("ForwardBarrel_dEvsPos_%d_%d",i,j),Form("dE vs pos in strip %d for fordet %d",j,i),330,0,1.1,300,0.,1.e3);
				// And position dependent
				for (Int_t k=0;k<16;k++) ForwardBarrel_dEforPos[i][j][k] = new TH1F(Form("ForwardBarrel_dEforPos_%d_%d_%d",i,j,k),Form("dE for pos %d in strip %d in fordet %d",k,j,i),500,0,5.);
	
				BackwardBarrel_dE[i][j] = new TH1F(Form("BackwardBarrel_dE_%d_%d",i,j),Form("dE for strip %d in backdet %d",j,i),300,0,15.);
				BackwardBarrel_dEforPos[i][j] = new TH1F(Form("BackwardBarrel_dEforPos_%d_%d",i,j),Form("dE for pos %d in backdet %d",j,i),500,0,15.);
	
				// E backdetector
				ForwardBarrel_dEvsE[i][j] = new TH2F(Form("ForwardBarrel_dEvsE_%d_%d",i,j),Form("dE vs Eback for strip %d in fordet %d",j,i),500,-1.,24.,500,0.,25.);
				ForwardBarrel_dEvsEcor[i][j] = new TH2F(Form("ForwardBarrel_dEvsEcor_%d_%d",i,j),Form("dE vs Eback for strip %d in fordet %d",j,i),500,-1.,24.,500,0.,25.);
				BackwardBarrel_dEvsE[i][j] = new TH2F(Form("BackwardBarrel_dEvsE_%d_%d",i,j),Form("dE vs Eback for strip %d in backdet %d",j,i),300,0.,6.,500,0.,10.);
	
				// Efor
				ForwardBarrel_EforvsPos[i][j] = new TH2F(Form("ForwardBarrel_EforvsPos_%d_%d",i,j),Form("Efor vs pos in strip %d for fordet %d",j,i),220,0,1.1,500,0,10.);
			}
	
			// Angle dependent histograms
			for (Int_t j =0;j<20;j++){
				// dEfile
				ForwardBarrel_dEangle[i][j] = new TH1F(Form("ForwardBarrel_dEangle_%d_%d",i,j),Form("dE around angle %d in fordet %d",j,i),500,0,1.);
				ForwardBarrel_dEvsEangle[i][j] = new TH2F(Form("ForwardBarrel_dEvsEangle_%d_%d",i,j),Form("dE vs E around angle %d in fordet %d",j,i),500,-1.,5.,500,0,25.);
				BackwardBarrel_dEvsEangle[i][j] = new TH2F(Form("BackwardBarrel_dEvsEangle_%d_%d",i,j),Form("dE vs E around angle %d in backdet %d",j,i),300,0,6.,500,0,10.);
			}
		}
	
		// Angle dependent histograms
		for (Int_t i=0;i<20;i++){
			Proton_Eex_Angle[i] = new TH1F(Form("Proton_Eex_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 50.,(i+1)*3. + 50.),500,-2.5,7.5);
			Proton_Eex_511_Angle[i] = new TH1F(Form("Proton_Eex_511_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 50.,(i+1)*3. + 50.),150,1.,4.);
			Proton_Eex_511_Angle_Background[i] = new TH1F(Form("Proton_Eex_511_Angle_Background_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 50.,(i+1)*3. + 50.),150,1.,4.);
			Proton_Eex_2100_Angle[i] = new TH1F(Form("Proton_Eex_2100_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 50.,(i+1)*3. + 50.),150,1.,4.);
			Proton_Eex_2100_Angle_Background[i] = new TH1F(Form("Proton_Eex_2100_Angle_Background_%d",i),Form("Background excitation for angles between %f and %f",i*3. + 50.,(i+1)*3. + 50.),150,1.,4.);
			Proton_Eex_2700_Angle[i] = new TH1F(Form("Proton_Eex_2700_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 50.,(i+1)*3. + 50.),150,1.,4.);
			Deuteron_Eex_Angle[i] = new TH1F(Form("Deuteron_Eex_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 60.,(i+1)*3. + 60.),500,-2.5,7.5);
			Deuteron_Eex_320_Angle[i] = new TH1F(Form("Deuteron_Eex_320_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 60.,(i+1)*3. + 60.),150,-1.5,1.5);
			Deuteron_Eex_320_Angle_Background[i] = new TH1F(Form("Deuteron_Eex_320_Angle_Background_%d",i),Form("Background excitaion for angles between %f and %f",i*3. + 60.,(i+1)*3. + 60.),150,-1.5,1.5);
			Triton_Eex_Angle[i] = new TH1F(Form("Triton_Eex_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 60.,(i+1)*3. + 60.),500,-2.5,7.5);
			Triton_Eex_3300_Angle[i] = new TH1F(Form("Triton_Eex_3300_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 60.,(i+1)*3. + 60.),150,2.,7.);
			Triton_Eex_3300_comp_Angle[i] = new TH1F(Form("Triton_Eex_3300_comp_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 60.,(i+1)*3. + 60.),150,2.,7.);
			Triton_Eex_2900_Angle[i] = new TH1F(Form("Triton_Eex_2900_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 60.,(i+1)*3. + 60.),150,2.,7.);
			Triton_Eex_2600_Angle[i] = new TH1F(Form("Triton_Eex_2600_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 60.,(i+1)*3. + 60.),150,2.,7.);
			Triton_Eex_6000_Angle[i] = new TH1F(Form("Triton_Eex_6000_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 60.,(i+1)*3. + 60.),150,4.,7.);
			Triton_Eex_6000_comp_Angle[i] = new TH1F(Form("Triton_Eex_6000_comp_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 60.,(i+1)*3. + 60.),150,4.,7.);
			NPT_Eex_Angle[i] = new TH1F(Form("NPT_Eex_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 80.,(i+1)*3. + 80.),500,-2.5,7.5);
			NPT_Eex_Angle_Back[i] = new TH1F(Form("NPT_Eex_Angle_Back_%d",i),Form("Excitaion energy for angles between %f and %f",i*2. + 140.,(i+1)*2. + 140.),150,-1.5,1.5);
			NPT_Eex_320_Angle[i] = new TH1F(Form("NPT_Eex_320_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 80.,(i+1)*3. + 80.),200,-2.,8.);
			NPT_Eex_320_Angle_Background[i] = new TH1F(Form("NPT_Eex_320_Angle_Background_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 80.,(i+1)*3. + 80.),200,-2.,8.);
			NPT_Eex_2100_Angle[i] = new TH1F(Form("NPT_Eex_2100_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 80.,(i+1)*3. + 80.),200,-2.,8.);
			NPT_Eex_2100_Angle_Background[i] = new TH1F(Form("NPT_Eex_2100_Angle_Background_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 80.,(i+1)*3. + 80.),200,-2.,8.);
			NPT_Eex_2600_Angle[i] = new TH1F(Form("NPT_Eex_2600_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 80.,(i+1)*3. + 80.),200,-2.,8.);
			NPT_Eex_2700_Angle[i] = new TH1F(Form("NPT_Eex_2700_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 80.,(i+1)*3. + 80.),200,-2.,8.);
			NPT_Eex_2900_Angle[i] = new TH1F(Form("NPT_Eex_2900_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 80.,(i+1)*3. + 80.),200,-2.,8.);
			NPT_Eex_3300_Angle[i] = new TH1F(Form("NPT_Eex_3300_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 80.,(i+1)*3. + 80.),200,-2.,8.);
			NPT_Eex_3300_comp_Angle[i] = new TH1F(Form("NPT_Eex_3300_comp_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 80.,(i+1)*3. + 80.),200,-2.,8.);
			NPT_Eex_3300_Angle_Background[i] = new TH1F(Form("NPT_Eex_3300_Angle_Background_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 80.,(i+1)*3. + 80.),200,-2.,8.);
			NPT_Eex_6000_comp_Angle[i] = new TH1F(Form("NPT_Eex_6000_comp_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 80.,(i+1)*3. + 80.),200,-2.,8.);
			NPT_Eex_6000_Angle_Background[i] = new TH1F(Form("NPT_Eex_6000_Angle_Background_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 80.,(i+1)*3. + 80.),200,-2.,8.);
			Beryllium_Eex_Angle_Stop_FB[i] = new TH1F(Form("Beryllium_Eex_Angle_Stop_FB_%d",i),Form("Excitaion energy for angles between %f and %f",i*2.5 + 25.,(i+1)*2.5 + 25.),500,-4.5,5.5);
			Beryllium_Eex_Angle_Stop_CD[i] = new TH1F(Form("Beryllium_Eex_Angle_Stop_CD_%d",i),Form("Excitaion energy for angles between %f and %f",i*1. + 8.,(i+1)*1. + 8.),500,-4.5,5.5);
		}
	
		// CD dependent histograms
		for(int i=0;i<4;i++){
			CD_dEtot[i] = new TH1F(Form("CD_dEtot_%d",i), Form("dE for CD %d",i),500,0,50.);
			CD_dEvsStrip[i] = new TH2F(Form("CD_dEvsStrip_%d",i), Form("deltaE vs Strip Nr of CD %d",i),16,0,16,250,0,50.);
			CD_dEvsRing[i] = new TH2F(Form("CD_dEvsRing_%d",i), Form("deltaE vs Ring Nr ofCD %d",i),16,0,16,250,0,50.);
			CD_StripvsRing[i] = new TH2F(Form("CD_StripvsRing_%d",i), Form("Strip vs Ring of CD %d",i),16,0,16,16,0,16);
			CD_EStripvsERing[i] = new TH2F(Form("CD_EStripvsERing_%d",i), Form("Strip vs Ring Energy of CD %d",i),250,0,50.,250,0,50.);
			CD_dEvsEtot[i] = new TH2F(Form("CD_dEvsEtot_%d",i), Form("dE vs E for CD %d",i),600,0.,30.,1000,0.,50.);
			CD_dEvsEcortot[i] = new TH2F(Form("CD_dEvsEcortot_%d",i), Form("Corrrected dE vs E for CD %d",i),600,0.,15.,1000,0.,25.);
			CD_Eback[i] = new TH1F(Form("CD_Eback_%d",i),Form("Eback for CD %d",i),600,0.,30.);
	
			// Angle
			CD_dEvsLabangle[i] = new TH2F(Form("CD_dEvsLabangle_%d",i),Form("dE vs Theta for CD %d",i),70,0,35,500,0.,50.);
	
	
	
			// Strip/ring dependent histograms
			for(int j=0;j<16;j++){
				CD_dEforStrip[i][j] = new TH1F(Form("CD_dEforStrip_%d_%d",i,j), Form("dE for strip %d in CD %d",j,i), 800,0,50);
				CD_dEforRing[i][j] = new TH1F(Form("CD_dEforRing_%d_%d",i,j), Form("dE for ring %d in CD %d",j,i), 800,0,50);
				CD_dEvsEStrip[i][j] = new TH2F(Form("CD_dEvsEStrip_%d_%d",i,j), Form("dE vs E for strip %d in CD %d",j,i),600,0.,30.,1000,0.,50.);
				CD_dEvsERing[i][j] = new TH2F(Form("CD_dEvsERing_%d_%d",i,j), Form("dE vs E for ring %d in CD %d",j,i),600,0.,30.,1000,0.,50.);
			}
		}
	
		// Gammas
		for (Int_t i=0;i<8;i++){
	
			// Crystals
			for (Int_t j=0;j<3;j++){
				Gamma_ECrystals1[i][j] = new TH1F(Form("Gamma_ECrystals_1_%d_%d",i,j),Form("Gamma energies in cluster %d and crystal %d",i,j),1000,0,4.);
				Gamma_ECrystals2[i][j] = new TH1F(Form("Gamma_ECrystals_2_%d_%d",i,j),Form("Gamma energies in cluster %d and crystal %d",i,j),1000,4.,8.);
	
				// Segments
				for (Int_t k=0;k<6;k++){
					Gamma_ESegments[i][j][k] = new TH1F(Form("Gamma_ESegments_%d_%d_%d",i,j,k),Form("Gamma energies in cluster %d,crystal %d and segment %d",i,j,k),800,0,8.);
					Proton_ESegments[i][j][k] = new TH1F(Form("Proton_ESegments_%d_%d_%d",i,j,k),Form("Gamma energies in cluster %d,crystal %d and segment %d",i,j,k),800,0,8.);
					Deuteron_ESegments[i][j][k] = new TH1F(Form("Deuteron_ESegments_%d_%d_%d",i,j,k),Form("Gamma energies in cluster %d,crystal %d and segment %d",i,j,k),800,0,8.);
	
				}
			}
		}
//	}
	// END HISTOGRAMS

	bool partfound;
	Double_t nentries = tr->GetEntries();
	Int_t nbytes = 0;
	Int_t status;

	Nd = 0;
	Np = 0;
	Nt = 0;
	Nd_FB = 0;
	Np_FB = 0;
	Nt_FB = 0;
	Nd_CD = 0;
	Np_CD = 0;
	Nt_CD = 0;
	tBe = 1;

	for(int i=0; i<nentries;i++){
		for (int idet=0;idet<3;idet++){
			for (int ie=0;ie<5;ie++){
				par[idet][ie] = 0;
				time[idet][ie] = 0;
				Eex[idet][ie] = 0;
				gamma[idet][ie] = 0;
				beta[idet][ie] = 0;
			}
		}
		status = tr->GetEvent(i);
		if(status == -1){
			cerr<<"Error occured, couldn't read entry "<<i<<" from tree "<<tr->GetName()<<" in file "<<tr->GetFile()->GetName()<<endl;
			return 5;
		}
		else if(status == 0){
			cerr<<"Error occured, entry "<<i<<" in tree "<<tr->GetName()<<" in file "<<tr->GetFile()->GetName()<<" doesn't exist"<<endl;
			return 6;
		}
		nbytes += status;


		partfound =false;

		// Forward barrel
		for(unsigned int j=0; j<FBarrel->size(); j++){
			if((*FBarrel)[j].GetMult() ==0) continue;

			if( (*FBarrel)[j].GetID() < 0 || (*FBarrel)[j].GetID() > 3 ){
				cerr<<"Error in entry "<<i<<": "<<j<<". forward barrel detector id is wrong: "<< (*FBarrel)[j].GetID()<<endl;
				continue;
			}

			if( (*FBarrel)[j].GetStripNr()[0] < 0 || (*FBarrel)[j].GetStripNr()[0] > 15 ){
				cerr<<"Error in entry "<<i<<": "<<(*FBarrel)[j].GetID()<<". forward barrel detector strip nr is wrong: "<<(*FBarrel)[j].GetStripNr()[0]<<endl;
				continue;
			}
			for (Int_t ih=0;ih<Histos.size();ih++){
				partfound=true;
				time[0][j] = (*FBarrel)[j].GetTime();
				Edet[0] = (*FBarrel)[j].GetEdet();
				
				// Cylindrical coordinates
				x = 50*((*FBarrel)[j].GetStripPos()[0]-.5) + xin*cos(pi/2.*(*FBarrel)[j].GetID()) + yin*sin(pi/2.*(*FBarrel)[j].GetID());
				y = 29. + xin*sin(pi/2.*(*FBarrel)[j].GetID()) - yin*cos(pi/2.*(*FBarrel)[j].GetID());
				z = 8. + 3.125*((*FBarrel)[j].GetStripNr()[0] + gRandom->Rndm()); 
				// Moving each detector
				for (Int_t ipad=0;ipad<4;ipad++){
					if ((*FBarrel)[j].GetID()==ipad){
						x += xpad[ipad];
						y += ypad[ipad];
						z += zpad[ipad];
					}
				}

				r = sqrt(pow(y,2) + pow(x,2));
				alpha = atan(x/y);
				// Spherical coordinates
				theta = atan(r/z);
				phi = alpha + pi/2.*((*FBarrel)[j].GetID()+1);
				if (phi< 0) phi += 2*pi;
				if (phi> 2*pi) phi -= 2*pi;
				// Detector angle
				thetadet = acos(sin(theta)*cos(alpha));
				
				// Energy loss in carbon foil for alphas
				if (strcmp(Histos[ih],"alpha")==0){
					Rafter = srim.RDL[4]->Eval((*FBarrel)[j].GetRear());
					Efor[0][4] = srim.EDL[4]->Eval(Rafter + .2/cos(thetadet));
					Rafter = srim.RMylar[4]->Eval(Efor[0][j]);
					Efor[0][4] = srim.EMylar[4]->Eval(Rafter + 11.57/cos(thetadet));
				}
				
				// Filling histos
				
				if (strcmp(Histos[ih],"sipos")==0){
					// Position histograms
					ForwardBarrel_Pos[(*FBarrel)[j].GetID()][(*FBarrel)[j].GetStripNr()[0]]->Fill((*FBarrel)[j].GetStripPos()[0]);
					ForwardBarrel_PosvsStrip[(*FBarrel)[j].GetID()]->Fill((*FBarrel)[j].GetStripNr()[0],(*FBarrel)[j].GetStripPos()[0]);
		
					// Angle histograms
					ForwardBarrel_Theta->Fill(theta/pi*180.);
					ForwardBarrel_Phi->Fill(phi/pi*180.);
					ForwardBarrel_Theta_det->Fill(thetadet);
		
					Hitpattern_ForwardBarreltotal->Fill(63*tan(theta)*cos(phi),63*tan(theta)*sin(phi));
				}

				if (strcmp(Histos[ih],"sie")==0){
					// dE histograms
					ForwardBarrel_dE[(*FBarrel)[j].GetID()][(*FBarrel)[j].GetStripNr()[0]]->Fill((*FBarrel)[j].GetRear()/1.e3);
					ForwardBarrel_dEvsPos[(*FBarrel)[j].GetID()][(*FBarrel)[j].GetStripNr()[0]]->Fill((*FBarrel)[j].GetStripPos()[0],(*FBarrel)[j].GetRear());
					ForwardBarrel_dEtot[(*FBarrel)[j].GetID()]->Fill((*FBarrel)[j].GetRear()/1.e3);
					ForwardBarrel_dEvsStrip[(*FBarrel)[j].GetID()]->Fill((*FBarrel)[j].GetStripNr()[0],(*FBarrel)[j].GetRear()/1.e3);
					ForwardBarrel_dEvsDetangle[(*FBarrel)[j].GetID()]->Fill(thetadet,(*FBarrel)[j].GetRear()/1.e3);
					for (Int_t M=0;M<16;M++){
						if (((*FBarrel)[j].GetStripPos()[0] > M*.06667)&&((*FBarrel)[j].GetStripPos()[0] < (M+1)*.06667)){
							ForwardBarrel_dEforPos[(*FBarrel)[j].GetID()][(*FBarrel)[j].GetStripNr()[0]][M]->Fill((*FBarrel)[j].GetRear()/1.e3);
						}
					}
					// dE for angle histograms
					for (Int_t k=0;k<20;k++){
						if ((thetadet > (k*.03 + .30)) && (thetadet < ((k+1)*.03 + .30))){
							ForwardBarrel_dEangle[(*FBarrel)[j].GetID()][k]->Fill((*FBarrel)[j].GetRear()/1.e3);
						}
						if ((thetadet > (k*.1 + .4)) && (thetadet < ((k+1)*.1 + .4))){
							ForwardBarrel_dEvsEangle[(*FBarrel)[j].GetID()][k]->Fill((*FBarrel)[j].GetEdet()/1.e3,(*FBarrel)[j].GetRear()/1.e3);
						}
					}
		
					// E back histograms
					ForwardBarrel_Eback[(*FBarrel)[j].GetID()]->Fill(Edet[0]/1.e3);
					ForwardBarrel_dEvsEtotal->Fill((*FBarrel)[j].GetEdet()/1.e3,(*FBarrel)[j].GetRear()/1.e3);
					ForwardBarrel_dEvsEtot[(*FBarrel)[j].GetID()]->Fill((*FBarrel)[j].GetEdet()/1.e3,(*FBarrel)[j].GetRear()/1.e3);
					ForwardBarrel_dEvsE[(*FBarrel)[j].GetID()][(*FBarrel)[j].GetStripNr()[0]]->Fill((*FBarrel)[j].GetEdet()/1.e3,(*FBarrel)[j].GetRear()/1.e3);
				}

				if (strcmp(Histos[ih],"alpha")==0){ 
					// Efor histograms
					ForwardBarrel_Efor[(*FBarrel)[j].GetID()]->Fill(Efor[0][4]/1.e3);
					ForwardBarrel_EforvsDetangle[(*FBarrel)[j].GetID()]->Fill(thetadet,Efor[0][4]/1.e3);
					ForwardBarrel_EforvsPos[(*FBarrel)[j].GetID()][(*FBarrel)[j].GetStripNr()[0]]->Fill((*FBarrel)[j].GetStripPos()[0],Efor[0][4]/1.e3);
				}
				
				// No punchthrough
				PT = srim.EDL[Apart]->Eval(140./cos(thetadet));
				if (((*FBarrel)[j].GetEdet() < 400.) && ((*FBarrel)[j].GetRear() > 500.) && ((*FBarrel)[j].GetRear()<PT)) {
//				if (((*FBarrel)[j].GetEdet() < 400.) && ((*FBarrel)[j].GetRear() > 500.)) {
					
					par[0][j] = 4;
					ID[0][j] = (*FBarrel)[j].GetID();
					thetap[0][j] = theta;
					phip[0][j] = phi;
					
					// ASSUMING ONE PARTICLE TYPE!!!!
					
					// Energy loss in carbon foil
					Rafter = srim.RMylar[Apart]->Eval((*FBarrel)[j].GetRear());
					Efor[0][j] = srim.EMylar[Apart]->Eval(Rafter + 11.57/cos(thetadet));
					if (Atarget==107){
						Rafter = srim.RT[20]->Eval(Efor[0][j]);
						Efor[0][j] = srim.ET[20]->Eval(Rafter + ttarget/(2.*cos(theta)));
					}
					else {
						Rafter = srim.RT[Apart]->Eval(Efor[0][j]);
						Efor[0][j] = srim.ET[Apart]->Eval(Rafter + ttarget/(2.*cos(theta)));
					}
					Px = sqrt(2*m[Apart]*Efor[0][j])*sin(theta)*cos(phi);
					Py = sqrt(2*m[Apart]*Efor[0][j])*sin(theta)*sin(phi);
					Pz = sqrt(2*m[Apart]*Efor[0][j])*cos(theta);
					Q = m[11] + m[Atarget] - m[Apart] - m[11+Atarget-Apart];
					EN = (pow(Pxin-Px,2) + pow(Pyin-Py,2) + pow(Pzin-Pz,2))/(2*m[11+Atarget-Apart]);
					
					Eex[0][j] = Ein + Q - Efor[0][j] - EN;
					// CM Angles
					vcm = (m[11]/(m[11]+m[Atarget]))*sqrt(2*Ein/m[11]);
					vout[0] = sqrt(2*Efor[0][j]/m[Apart])*cos(theta);
					vout[1] = sqrt(2*Efor[0][j]/m[Apart])*sin(theta);
					theta_cm[0][j] = atan(vout[1]/(vout[0]-vcm));
					if (theta_cm[0][j] < 0) theta_cm[0][j] += pi;

					// Relativistic factors
					gamma[0][j] = (Ein + (m[11] + m[Atarget] - m[11+Atarget-Apart] - m[Apart]) - Eex[0][j] - Efor[0][j] + m[11+Atarget-Apart])/m[11+Atarget-Apart];
					beta[0][j] = sqrt(1. - 1./(gamma[0][j]*gamma[0][j]));
					
					
					if (strcmp(Histos[ih],"part")==0){ 
					
					
						if (((*FBarrel)[j].GetID()==2)||((*FBarrel)[j].GetID()==3)){
							ForwardBarrel_dEvsLabangle[(*FBarrel)[j].GetID()]->Fill(theta/pi*180.,Efor[0][j]/1.e3);
							ForwardBarrel_dEvsLabangletotal->Fill(theta/pi*180.,Efor[0][j]/1.e3);
							ForwardBarrel_dEvsCmangletotal->Fill(theta_cm[0][j]/pi*180.,Efor[0][j]/1.e3);
							All_dEvsLabangletotal->Fill(theta,(*FBarrel)[j].GetRear());
							NPT_Eextotal->Fill(Eex[0][j]/1.e3);
							for (Int_t iAngle=0; iAngle<20.;iAngle++)
							{
								if ((theta_cm[0][j]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[0][j]<((iAngle+1)*3./180.*pi + 80./180.*pi))){
									NPT_Eex_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
										if (sim){
											NPT_Eex_320_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
											NPT_Eex_2100_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
											NPT_Eex_2600_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
											NPT_Eex_2700_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
											NPT_Eex_2900_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
											NPT_Eex_3300_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
											NPT_Eex_6000_comp_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
										}
								}
								if ((theta>(iAngle*2.5/180.*pi + 25./180.*pi)) && (theta<((iAngle+1)*2.5/180.*pi + 25./180.*pi)) && (Apart==11)){
									Beryllium_Eex_Angle_Stop_FB[iAngle]->Fill(Eex[0][j]/1.e3);
								}
							
							}
//							printf("%f, %f \n",(*FBarrel)[j].GetRear(),Be_cut[1]->Eval(theta));
							if (((*FBarrel)[j].GetRear() > Be_cut[0]->Eval(theta)) && ((*FBarrel)[j].GetRear() < Be_cut[1]->Eval(theta))){
								Be11_Angle->Fill(theta/pi*180.);
								Be11_pr_Time->Fill(EbisTime*25.e-9);
							}
						}
					}
				}
	
				// Particle identification
				if (((*FBarrel)[j].GetEdet() > 400.) && ((*FBarrel)[j].GetRear() > 500.)){
					Ecor = (*FBarrel)[j].GetEdet() + (1. - cos(thetadet))*(*FBarrel)[j].GetRear();
					dEcor = (*FBarrel)[j].GetRear()*cos(thetadet);
					
					if ((strcmp(Histos[ih],"part")==0)||(strcmp(Histos[ih],"sie")==0)){ 
						ForwardBarrel_dEvsEcortotal->Fill(Ecor/1.e3,dEcor/1.e3);
						ForwardBarrel_dEvsEcortot[(*FBarrel)[j].GetID()]->Fill(Ecor/1.e3,dEcor/1.e3);
						ForwardBarrel_dEvsEcor[(*FBarrel)[j].GetID()][(*FBarrel)[j].GetStripNr()[0]]->Fill(Ecor/1.e3,dEcor/1.e3);
					}
	
					// Protons
					if ((dEcor > For_cut[0][(*FBarrel)[j].GetStripNr()[0]]->Eval(Ecor)) && (dEcor < For_cut[1][(*FBarrel)[j].GetStripNr()[0]]->Eval(Ecor))){
						par[0][j] = 1;
						ID[0][j] = (*FBarrel)[j].GetID();
						Np++;
						Np_FB++;
						thetap[0][j] = theta;
						phip[0][j] = phi;
						// Energy loss in carbon foil
						Rafter = srim.RMylar[1]->Eval((*FBarrel)[j].GetEdet() + (*FBarrel)[j].GetRear());
						Efor[0][j] = srim.EMylar[1]->Eval(Rafter + 11.57/cos(thetadet));
						Rafter = srim.RT[1]->Eval(Efor[0][j]);
						Efor[0][j] = srim.ET[1]->Eval(Rafter + ttarget/(2.*cos(theta)));
						Px = sqrt(2*m[1]*Efor[0][j])*sin(theta)*cos(phi);
						Py = sqrt(2*m[1]*Efor[0][j])*sin(theta)*sin(phi);
						Pz = sqrt(2*m[1]*Efor[0][j])*cos(theta);
						Q = m[11] + m[Atarget] - m[1] - m[11+Atarget-1];
						EN = (pow(Pxin-Px,2) + pow(Pyin-Py,2) + pow(Pzin-Pz,2))/(2*m[11+Atarget-1]);
						
						Eex[0][j] = Ein + Q - Efor[0][j] - EN;
						// CM Angles
						vcm = (m[11]/(m[11]+m[Atarget]))*sqrt(2*Ein/m[11]);
						vout[0] = sqrt(2*Efor[0][j]/m[1])*cos(theta);
						vout[1] = sqrt(2*Efor[0][j]/m[1])*sin(theta);
						theta_cm[0][j] = atan(vout[1]/(vout[0]-vcm));
						if (theta_cm[0][j] < 0) theta_cm[0][j] += pi;
						// Relativistic factors
						gamma[0][j] = (Ein + (m[11] + m[Atarget] - m[11+Atarget-1] - m[1]) - Eex[0][j] - Efor[0][j] + m[11+Atarget-1])/m[11+Atarget-1];
						//gamma[0][j] = (Ein + (m[22] + m[2] - m[23] - m[1]) - Eex[0][j] - Efor[0][j] + m[23])/m[23];
						beta[0][j] = sqrt(1. - 1./(gamma[0][j]*gamma[0][j]));
				
						if (strcmp(Histos[ih],"part")==0){	
							// Filling histograms
							ForwardBarrel_dEvsEcorp->Fill(Ecor/1.e3,dEcor/1.e3);
							Proton_EtotvsLabangle[(*FBarrel)[j].GetID()]->Fill(theta/pi*180.,Efor[0][j]/1.e3);
							Proton_EtotvsLabangletotal->Fill(theta/pi*180.,Efor[0][j]/1.e3);
							Proton_EtotvsCmangletotal->Fill(theta_cm[0][j]/pi*180.,Efor[0][j]/1.e3);
							Proton_Eex_ForwardBarrel[(*FBarrel)[j].GetID()]->Fill(Eex[0][j]/1.e3);
							if ((*FBarrel)[j].GetID()!=0) Proton_Eextotal->Fill(Eex[0][j]/1.e3);
							Proton_EexvsLabangle_ForwardBarrel[(*FBarrel)[j].GetID()]->Fill(theta/pi*180.,Eex[0][j]/1.e3);
							
							if ((*FBarrel)[j].GetID()!=0){
								for (Int_t iAngle=0; iAngle<20.;iAngle++)
								{
									if ((theta_cm[0][j]>(iAngle*3./180.*pi + 50./180.*pi)) && (theta_cm[0][j]<((iAngle+1)*3./180.*pi + 50./180.*pi))){ 
										Proton_Eex_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
										if (sim){
											Proton_Eex_511_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
											Proton_Eex_2100_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
											Proton_Eex_2700_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
										}
									}
								}
							}
						}
					}
	
					// Deuterons
					if ((dEcor > For_cut[1][(*FBarrel)[j].GetStripNr()[0]]->Eval(Ecor)) && (dEcor < For_cut[2][(*FBarrel)[j].GetStripNr()[0]]->Eval(Ecor))){
						par[0][j] = 2;
						ID[0][j] = (*FBarrel)[j].GetID();
						Nd++;
						Nd_FB++;
						Deuteron_pr_Time->Fill(EbisTime*25.e-9);
					
						thetap[0][j] = theta;
						phip[0][j] = phi;
						// Energy loss in carbon foil
						Rafter = srim.RMylar[2]->Eval((*FBarrel)[j].GetEdet() + (*FBarrel)[j].GetRear());
						Efor[0][j] = srim.EMylar[2]->Eval(Rafter + 11.57/cos(thetadet));
						Rafter = srim.RT[2]->Eval(Efor[0][j]);
						Efor[0][j] = srim.ET[2]->Eval(Rafter + ttarget/(2.*cos(theta)));
						Px = sqrt(2*m[2]*Efor[0][j])*sin(theta)*cos(phi);
						Py = sqrt(2*m[2]*Efor[0][j])*sin(theta)*sin(phi);
						Pz = sqrt(2*m[2]*Efor[0][j])*cos(theta);
						Q = m[11] + m[Atarget] - m[2] - m[11+Atarget-2];
						EN = (pow(Pxin-Px,2) + pow(Pyin-Py,2) + pow(Pzin-Pz,2))/(2*m[11+Atarget-2]);
						Eex[0][j] = Ein + Q - Efor[0][j] - EN;
						// CM Angles
						vcm = (m[11]/(m[11]+m[Atarget]))*sqrt(2*Ein/m[11]);
						vout[0] = sqrt(2*Efor[0][j]/m[2])*cos(theta);
						vout[1] = sqrt(2*Efor[0][j]/m[2])*sin(theta);
						theta_cm[0][j] = atan(vout[1]/(vout[0]-vcm));
						if (theta_cm[0][j] < 0) theta_cm[0][j] += pi;
						// Relativistic factors
						gamma[0][j] = (Ein + (m[11] + m[Atarget] - m[11+Atarget-2] - m[2]) - Eex[0][j] - Efor[0][j] + m[11+Atarget-2])/m[11+Atarget-2];
						//gamma[0][j] = (Ein + (m[22] + m[2] - m[22] - m[2]) - Eex[0][j] - Efor[0][j] + m[22])/m[22];
						beta[0][j] = sqrt(1. - 1./(gamma[0][j]*gamma[0][j]));
					
						if (strcmp(Histos[ih],"part")==0){
							// Filling histograms
							ForwardBarrel_dEvsEcord->Fill(Ecor/1.e3,dEcor/1.e3);
							Deuteron_EtotvsLabangle[(*FBarrel)[j].GetID()]->Fill(theta/pi*180.,Efor[0][j]/1.e3);
							Deuteron_EtotvsLabangletotal->Fill(theta/pi*180.,Efor[0][j]/1.e3);
							Deuteron_EtotvsCmangletotal->Fill(theta_cm[0][j]/pi*180.,Efor[0][j]/1.e3);
							Deuteron_Eex_ForwardBarrel[(*FBarrel)[j].GetID()]->Fill(Eex[0][j]/1.e3);
							if ((*FBarrel)[j].GetID()!=0) Deuteron_Eextotal->Fill(Eex[0][j]/1.e3);
							Deuteron_EexvsLabangle_ForwardBarrel[(*FBarrel)[j].GetID()]->Fill(theta/pi*180.,Eex[0][j]/1.e3);
		
							if ((*FBarrel)[j].GetID()!=0){
								for (Int_t iAngle=0; iAngle<20.;iAngle++)
								{
									if ((theta_cm[0][j]>(iAngle*3./180.*pi + 60./180.*pi)) && (theta_cm[0][j]<((iAngle+1)*3./180.*pi + 60./180.*pi))){
										Deuteron_Eex_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
										if (sim){
										Deuteron_Eex_320_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
										}
									}
								}
							}
						}
					}
	
					// Tritons
					if ((dEcor > For_cut[2][(*FBarrel)[j].GetStripNr()[0]]->Eval(Ecor)) && (dEcor < For_cut[3][(*FBarrel)[j].GetStripNr()[0]]->Eval(Ecor))){
						par[0][j] = 3;
						ID[0][j] = (*FBarrel)[j].GetID();
						Nt++;
						Nt_FB++;
						thetap[0][j] = theta;
						phip[0][j] = phi;
						// Energy loss in carbon foil
						Rafter = srim.RMylar[3]->Eval((*FBarrel)[j].GetEdet() + (*FBarrel)[j].GetRear());
						Efor[0][j] = srim.EMylar[3]->Eval(Rafter + 11.57/cos(thetadet));
						Rafter = srim.RT[3]->Eval(Efor[0][j]);
						Efor[0][j] = srim.ET[3]->Eval(Rafter + ttarget/(2.*cos(theta)));
						Px = sqrt(2*m[3]*Efor[0][j])*sin(theta)*cos(phi);
						Py = sqrt(2*m[3]*Efor[0][j])*sin(theta)*sin(phi);
						Pz = sqrt(2*m[3]*Efor[0][j])*cos(theta);
						Q = m[11] + m[Atarget] - m[3] - m[11+Atarget-3];
						EN = (pow(Pxin-Px,2) + pow(Pyin-Py,2) + pow(Pzin-Pz,2))/(2*m[11+Atarget-3]);
						
						Eex[0][j] = Ein + Q - Efor[0][j] - EN;
//						printf("%f, %f \n",Eex[0][j],Q);
						// CM Angles
						vcm = (m[11]/(m[11]+m[Atarget]))*sqrt(2*Ein/m[11]);
						vout[0] = sqrt(2*Efor[0][j]/m[3])*cos(theta);
						vout[1] = sqrt(2*Efor[0][j]/m[3])*sin(theta);
						theta_cm[0][j] = atan(vout[1]/(vout[0]-vcm));
						if (theta_cm[0][j] < 0) theta_cm[0][j] += pi;
						// Relativistic factors
						gamma[0][j] = (Ein + (m[11] + m[Atarget] - m[11+Atarget-3] - m[3]) - Eex[0][j] - Efor[0][j] + m[11+Atarget-3])/m[11+Atarget-3];
						beta[0][j] = sqrt(1. - 1./(gamma[0][j]*gamma[0][j]));
							
						
						if (strcmp(Histos[ih],"part")==0){
							// Filling histograms
							ForwardBarrel_dEvsEcort->Fill(Ecor/1.e3,dEcor/1.e3);
							Triton_EtotvsLabangle[(*FBarrel)[j].GetID()]->Fill(theta/pi*180.,Efor[0][j]/1.e3);
							Triton_EtotvsLabangletotal->Fill(theta/pi*180.,Efor[0][j]/1.e3);
							Triton_EtotvsCmangletotal->Fill(theta_cm[0][j]/pi*180.,Efor[0][j]/1.e3);
							Triton_Eex_ForwardBarrel[(*FBarrel)[j].GetID()]->Fill(Eex[0][j]/1.e3);
							if ((*FBarrel)[j].GetID()!=0) Triton_Eextotal->Fill(Eex[0][j]/1.e3);
							Triton_EexvsLabangle_ForwardBarrel[(*FBarrel)[j].GetID()]->Fill(theta/pi*180.,Eex[0][j]/1.e3);
		
							if ((*FBarrel)[j].GetID()!=0){
								for (Int_t iAngle=0; iAngle<20.;iAngle++)
								{
									if ((theta_cm[0][j]>(iAngle*3./180.*pi + 60./180.*pi)) && (theta_cm[0][j]<((iAngle+1)*3./180.*pi + 60./180.*pi))){
										Triton_Eex_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
										if (sim){
											Triton_Eex_3300_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
											Triton_Eex_3300_comp_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
											Triton_Eex_2600_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
											Triton_Eex_2900_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
											Triton_Eex_6000_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
											Triton_Eex_6000_comp_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
										}
									}
								}
							}
						}	
					}
				}
			}
		}
		// End forward barrel

		// Backward barrel
		for(unsigned int j=0; j<BBarrel->size(); j++){
			if((*BBarrel)[j].GetMult() ==0)
				continue;
			if( (*BBarrel)[j].GetID() < 0 || (*BBarrel)[j].GetID() > 3 ){
				cerr<<"Error in entry "<<i<<": "<<j<<". backward barrel detector id is wrong: "<<(*BBarrel)[j].GetID()<<endl;
				continue;
			}
			if( (*BBarrel)[j].GetStripNr()[0] < 0 || (*BBarrel)[j].GetStripNr()[0] > 15 ){
				cerr<<"Error in entry "<<i<<": "<<j<<". backward barrel detector strip nr is wrong: "<<(*BBarrel)[j].GetStripNr()[0]<<endl;
				continue;
			}
			for (Int_t ih=0;ih<Histos.size();ih++){
				time[1][j] = (*BBarrel)[j].GetTime();
				Edet[1] = (*BBarrel)[j].GetEdet();
	
				// Cylindrical coordinates
				x = 50*((*BBarrel)[j].GetStripPos()[0]-.5) + xin*cos(pi/2.*(*BBarrel)[j].GetID()) + yin*sin(pi/2.*(*BBarrel)[j].GetID());
				y = 29. + xin*sin(pi/2.*(*BBarrel)[j].GetID()) - yin*cos(pi/2.*(*BBarrel)[j].GetID());
				z = 8. + 3.125*((*BBarrel)[j].GetStripNr()[0] + gRandom->Rndm()); 
				r = sqrt(pow(y,2) + pow(x,2));
				alpha = atan(x/y);
				
				// Spherical coordinates
				theta = atan(r/z);
				phi = alpha + pi/2.*((*BBarrel)[j].GetID()+1);
				if (phi< 0) phi += 2*pi;
				if (phi> 2*pi) phi -= 2*pi;
				// Detector angle
				thetadet = acos(sin(theta)*cos(alpha));
				
				// Filling histos
				if (strcmp(Histos[ih],"sipos")==0){
					// Position histograms
					BackwardBarrel_Pos[(*BBarrel)[j].GetID()][(*BBarrel)[j].GetStripNr()[0]]->Fill((*BBarrel)[j].GetStripPos()[0]);
					BackwardBarrel_PosvsStrip[(*BBarrel)[j].GetID()]->Fill((*BBarrel)[j].GetStripNr()[0],(*BBarrel)[j].GetStripPos()[0]);
					// Angle histograms
					BackwardBarrel_Theta->Fill(theta/pi*180.);
					BackwardBarrel_Phi->Fill(phi/pi*180.);
				}

				if (strcmp(Histos[ih],"sie")==0){
					// dE histograms
					BackwardBarrel_dE[(*BBarrel)[j].GetID()][(*BBarrel)[j].GetStripNr()[0]]->Fill((*BBarrel)[j].GetRear()/1.e3);
					BackwardBarrel_dEvsPos[(*BBarrel)[j].GetID()]->Fill((*BBarrel)[j].GetStripPos()[0],(*BBarrel)[j].GetRear());
					BackwardBarrel_dEtot[(*BBarrel)[j].GetID()]->Fill((*BBarrel)[j].GetRear()/1.e3);
					BackwardBarrel_dEvsStrip[(*BBarrel)[j].GetID()]->Fill((*BBarrel)[j].GetStripNr()[0],(*BBarrel)[j].GetRear()/1.e3);
					for (Int_t M=0;M<16;M++){
						if (((*BBarrel)[j].GetStripPos()[0] > M*.06667)&&((*BBarrel)[j].GetStripPos()[0] < (M+1)*.06667)){
							BackwardBarrel_dEforPos[(*BBarrel)[j].GetID()][M]->Fill((*BBarrel)[j].GetRear()/1.e3);
						}
					}
		
					// E back histograms
					BackwardBarrel_Eback[(*BBarrel)[j].GetID()]->Fill(Edet[1]/1.e3);
					BackwardBarrel_dEvsEtotal->Fill((*BBarrel)[j].GetEdet()/1.e3,(*BBarrel)[j].GetRear()/1.e3);
					BackwardBarrel_dEvsEtot[(*BBarrel)[j].GetID()]->Fill((*BBarrel)[j].GetEdet()/1.e3,(*BBarrel)[j].GetRear()/1.e3);
					BackwardBarrel_dEvsE[(*BBarrel)[j].GetID()][(*BBarrel)[j].GetStripNr()[0]]->Fill((*BBarrel)[j].GetEdet()/1.e3,(*BBarrel)[j].GetRear()/1.e3);
					// dE for angle histograms
					for (Int_t k=0;k<6;k++){
						if ((thetadet > (k*.1 + .4)) && (thetadet < ((k+1)*.1 + .4))){
							BackwardBarrel_dEvsEangle[(*BBarrel)[j].GetID()][k]->Fill((*BBarrel)[j].GetEdet()/1.e3,(*BBarrel)[j].GetRear()/1.e3);
						}
					}
				}
				
				// No punchthrough
				if (((*BBarrel)[j].GetEdet() < 200.) && ((*BBarrel)[j].GetRear() > 50.)){
					par[1][j] = 4;
					
					// Energy loss in carbon foil
					Rafter = srim.RT[1]->Eval((*BBarrel)[j].GetRear());
//					Efor[1][j] = srim.EDL[1]->Eval(Rafter + .5/cos(thetadet));
//					Rafter = srim.RT[1]->Eval(Efor[1][j]);
					Efor[1][j] = srim.ET[1]->Eval(Rafter + ttarget/(2.*cos(theta)));
					Px = sqrt(2*m[1]*Efor[1][j])*sin(pi-theta)*cos(phi);
					Py = sqrt(2*m[1]*Efor[1][j])*sin(pi-theta)*sin(phi);
					Pz = sqrt(2*m[1]*Efor[1][j])*cos(pi-theta);
					Q = m[11] + m[Atarget] - m[1] - m[11+Atarget-1];
					EN = (pow(Pxin-Px,2) + pow(Pyin-Py,2) + pow(Pzin-Pz,2))/(2*m[11+Atarget-1]);
					Eex[1][j] = Ein + Q - Efor[1][j] - EN;
					
					// CM Angles
					vcm = (m[11]/(m[11]+m[Atarget]))*sqrt(2*Ein/m[11]);
					vout[0] = sqrt(2*Efor[1][j]/m[1])*cos(pi-theta);
					vout[1] = sqrt(2*Efor[1][j]/m[1])*sin(pi-theta);
					theta_cm[1][j] = atan(vout[1]/(vout[0]-vcm));
					if (theta_cm[1][j] < 0) theta_cm[1][j] += pi;
					
					// Relativistic factors
					gamma[1][j] = (Ein + (m[11] + m[Atarget] - m[11+Atarget-1] - m[1]) - Eex[1][j] - Efor[1][j] + m[11+Atarget-1])/m[11+Atarget-1];
					beta[1][j] = sqrt(1. - 1./(gamma[1][j]*gamma[1][j]));
					
					if (strcmp(Histos[ih],"part")==0){
						BackwardBarrel_dEvsLabangle[(*BBarrel)[j].GetID()]->Fill(180. - theta/pi*180.,(*BBarrel)[j].GetRear()/1.e3);
						BackwardBarrel_dEvsLabangletotal->Fill(180. - theta/pi*180.,(*BBarrel)[j].GetRear()/1.e3);
						NPT_Eex_BackwardBarrel[(*BBarrel)[j].GetID()]->Fill(Eex[1][j]/1.e3);
						NPT_Eextotal_Back->Fill(Eex[1][j]/1.e3);
						
						if ((*BBarrel)[j].GetID()==3){
							for (Int_t iAngle=0; iAngle<20.;iAngle++)
							{
								if ((theta_cm[1][j]>(iAngle*2./180.*pi + 140./180.*pi)) && (theta_cm[1][j]<((iAngle+1)*2./180.*pi + 140./180.*pi))){
									NPT_Eex_Angle_Back[iAngle]->Fill(Eex[1][j]/1.e3);
								}
							}
						}
					}
				}
	
				// Particle identification
				if (((*BBarrel)[j].GetEdet() > 200.)){
					Ecor = (*BBarrel)[j].GetEdet() + (1. - cos(thetadet))*(*BBarrel)[j].GetRear();
					dEcor = (*BBarrel)[j].GetRear()*cos(thetadet);
					if ((strcmp(Histos[ih],"part")==0)||(strcmp(Histos[ih],"sie")==0)){
						BackwardBarrel_dEvsEcortot[(*BBarrel)[j].GetID()]->Fill(Ecor/1.e3,dEcor/1.e3);
						BackwardBarrel_dEvsEcortotal->Fill(Ecor/1.e3,dEcor/1.e3);
					}
					// Protons
					if ((dEcor > Back_cut[0]->Eval(Ecor)) && (dEcor < Back_cut[1]->Eval(Ecor))){
						par[1][j] = 1;
						// Energy loss
						Rafter = srim.RDL[1]->Eval((*BBarrel)[j].GetEdet() + (*BBarrel)[j].GetRear());
						Efor[1][j] = srim.EDL[1]->Eval(Rafter + .5/cos(thetadet));
						Rafter = srim.RT[1]->Eval(Efor[1][j]);
						Efor[1][j] = srim.ET[1]->Eval(Rafter + ttarget/(2.*cos(theta)));
						Px = sqrt(2*m[1]*Efor[1][j])*sin(pi-theta)*cos(phi);
						Py = sqrt(2*m[1]*Efor[1][j])*sin(pi-theta)*sin(phi);
						Pz = sqrt(2*m[1]*Efor[1][j])*cos(pi-theta);
						Q = m[11] + m[Atarget] - m[1] - m[11+Atarget-1];
						EN = (pow(Pxin-Px,2) + pow(Pyin-Py,2) + pow(Pzin-Pz,2))/(2*m[11+Atarget-1]);
						
						Eex[1][j] = Ein + Q - Efor[1][j] - EN;
						// CM Angles
						vcm = (m[11]/(m[11]+m[Atarget]))*sqrt(2*Ein/m[11]);
						vout[0] = sqrt(2*Efor[1][j]/m[1])*cos(theta);
						vout[1] = sqrt(2*Efor[1][j]/m[1])*sin(theta);
						theta_cm[1][j] = atan(vout[1]/(vout[0]-vcm));
						if (theta_cm[1][j] < 0) theta_cm[1][j] += pi;
						// Relativistic factors
						gamma[1][j] = (Ein + Q - Eex[1][j] - Efor[1][j] + m[11+Atarget-1])/m[11+Atarget-1];
						beta[1][j] = sqrt(1. - 1./(gamma[1][j]*gamma[1][j]));
	
						if (strcmp(Histos[ih],"part")==0){
							BackwardBarrel_dEvsEcorp->Fill(Ecor/1.e3,dEcor/1.e3);
							Proton_Back_EtotvsLabangletotal->Fill(-theta/pi*180. + 180.,Efor[1][j]/1.e3);
						}
					}
				}
			}
		}
		// End backward barrel

		// Forward CD
		for(unsigned int j=0; j<FCD->size(); j++){
			if( (*FCD)[j].GetID() < 0 || (*FCD)[j].GetID() >= 4 ){
				cerr<<"Error in entry "<<i<<": "<<j<<". backward barrel detector id is wrong: "<<(*FCD)[j].GetID()<<endl;
				continue;
			}

			time[2][j] = (*FCD)[j].GetTime();
			Edet[2] = (*FCD)[j].GetEdet();

			// Filling histos

			// Ring dependent histograms
			for( unsigned int k=0;k<(*FCD)[j].GetRingNr().size();k++){
				if(verbose) cout << "ring nr " << (*FCD)[j].GetRingNr()[k] << " with energy " << (*FCD)[j].GetRingEnergy()[k] << endl;
				for (Int_t ih=0;ih<Histos.size();ih++){
					
					if (strcmp(Histos[ih],"sie")==0){
						// dE
						CD_dEvsRing[(*FCD)[j].GetID()]->Fill( (*FCD)[j].GetRingNr()[k], (*FCD)[j].GetRingEnergy()[k]/1.e3);
						CD_dEforRing[(*FCD)[j].GetID()][(*FCD)[j].GetRingNr()[k]]->Fill((*FCD)[j].GetRingEnergy()[k]/1.e3);
		
						// E back histograms
						CD_Eback[(*FCD)[j].GetID()]->Fill((*FCD)[j].GetEdet()/1.e3);
						CD_dEvsEtot[(*FCD)[j].GetID()]->Fill((*FCD)[j].GetEdet()/1.e3,(*FCD)[j].GetRingEnergy()[k]/1.e3);
						CD_dEvsERing[(*FCD)[j].GetID()][(*FCD)[j].GetRingNr()[k]]->Fill((*FCD)[j].GetEdet()/1.e3,(*FCD)[j].GetRingEnergy()[k]/1.e3);
					}
					
					// Strip and ring dependent histograms
					for(unsigned int l=0;l<(*FCD)[j].GetStripNr().size();l++){
						if(verbose) cout << "strip nr " << (*FCD)[j].GetStripNr()[l] << " with energy " << (*FCD)[j].GetStripEnergy()[l] << endl;
	
						if (k==0){
							if (strcmp(Histos[ih],"sie")==0){
								CD_dEvsStrip[(*FCD)[j].GetID()]->Fill( (*FCD)[j].GetStripNr()[l], (*FCD)[j].GetStripEnergy()[l]/1.e3);
								CD_dEforStrip[(*FCD)[j].GetID()][(*FCD)[j].GetStripNr()[l]]->Fill((*FCD)[j].GetStripEnergy()[l]/1.e3);
							}
						}
							
						// Angle
						// Values from centrum
						r = 9. + (gRandom->Rndm() + (*FCD)[j].GetRingNr()[k])*2.;	
						z = 63.;
						theta = atan(r/z);

						if (((*FCD)[j].GetStripNr()[l]>=4)&&((*FCD)[j].GetStripNr()[l]<12)){
							phi = .0593411*(4+ 2*((*FCD)[j].GetStripNr()[l]-4 + gRandom->Rndm())) - .71209;
						}
						else if ((*FCD)[j].GetStripNr()[l]<4){
							phi = .0593411*((*FCD)[j].GetStripNr()[l] + gRandom->Rndm()) - .71209;
						}
						else{
							phi = .0593411*((*FCD)[j].GetStripNr()[l]+8 + gRandom->Rndm()) - .71209;
						}
						
						// Moving to beamspot
						x = r*sin(phi) + xin*cos(pi/2.*(*FCD)[j].GetID()) + yin*sin(pi/2.*(*FCD)[j].GetID()); 
						y = r*cos(phi) + xin*sin(pi/2.*(*FCD)[j].GetID()) - yin*cos(pi/2.*(*FCD)[j].GetID()); 
						r = sqrt(pow(y,2)+pow(x,2));
						theta = atan(r/z);
						phi = atan(x/y);
						phi = phi + pi/2*((*FCD)[j].GetID()+1);
						if (phi>2*pi) phi -= 2*pi;
						if (phi<0) phi += 2*pi;
						
						//Fill histograms
						if (strcmp(Histos[ih],"sipos")==0){
							CD_Phi->Fill(phi/pi*180.);
							CD_Theta->Fill(theta/pi*180.);
						}
	
						if (TMath::Abs((*FCD)[j].GetRingEnergy()[k] - (*FCD)[j].GetStripEnergy()[l]) < 500){
							
							if (strcmp(Histos[ih],"sipos")==0) Hitpattern_CDtotal->Fill(63*tan(theta)*cos(phi),63*tan(theta)*sin(phi));
							
							if (strcmp(Histos[ih],"sipos")==0) CD_StripvsRing[(*FCD)[j].GetID()]->Fill( (*FCD)[j].GetStripNr()[l], (*FCD)[j].GetRingNr()[k] );
						
							if (strcmp(Histos[ih],"sie")==0){
								CD_EStripvsERing[(*FCD)[j].GetID()]->Fill( (*FCD)[j].GetStripEnergy()[l]/1.e3, (*FCD)[j].GetRingEnergy()[k]/1.e3);
								CD_dEtot[(*FCD)[j].GetID()]->Fill((*FCD)[j].GetStripEnergy()[l]/1.e3);
								CD_dEvsEtot[(*FCD)[j].GetID()]->Fill((*FCD)[j].GetEdet()/1.e3,(*FCD)[j].GetStripEnergy()[l]/1.e3);
								CD_dEvsEtotal->Fill((*FCD)[j].GetEdet()/1.e3,(*FCD)[j].GetStripEnergy()[l]/1.e3);
								CD_dEvsEStrip[(*FCD)[j].GetID()][(*FCD)[j].GetStripNr()[l]]->Fill((*FCD)[j].GetEdet()/1.e3,(*FCD)[j].GetStripEnergy()[l]/1.e3);
							}

							// No Punchthrough
							PT = srim.EDL[Apart]->Eval(500/cos(theta));
//							if ((*FCD)[j].GetEdet() < CDThreshold[(*FCD)[j].GetID()] && ((*FCD)[j].GetRingEnergy()[k]<PT)){
							if ((*FCD)[j].GetEdet() < CDThreshold[(*FCD)[j].GetID()]){
								par[2][j] = 4;
								thetap[2][j] = theta;
							// ASSUMING ONE PARTICLE TYPE!!!!
							if (Apart!=4){
							Rafter = srim.RDL[11+Atarget-Apart]->Eval((*FCD)[j].GetRingEnergy()[k]);
							Efor[2][j] = srim.EDL[11+Atarget-Apart]->Eval(Rafter + .7/cos(theta));
							if (Atarget==107){
								Rafter = srim.RT[20]->Eval(Efor[2][j]);
								Efor[2][j] = srim.ET[20]->Eval(Rafter + ttarget/(2.*cos(theta)));
							}
							else {
								Rafter = srim.RT[11+Atarget-Apart]->Eval(Efor[2][j]);
								Efor[2][j] = srim.ET[11+Atarget-Apart]->Eval(Rafter + ttarget/(2.*cos(theta)));
							}
							}
							Px = sqrt(2*m[Apart]*Efor[2][j])*sin(theta)*cos(phi);
							Py = sqrt(2*m[Apart]*Efor[2][j])*sin(theta)*sin(phi);
							Pz = sqrt(2*m[Apart]*Efor[2][j])*cos(theta);
							Q = m[11] + m[Atarget] - m[Apart] - m[11+Atarget-Apart];
							EN = (pow(Pxin-Px,2) + pow(Pyin-Py,2) + pow(Pzin-Pz,2))/(2*m[11+Atarget-Apart]);
							
							Eex[2][j] = Ein + Q - Efor[2][j] - EN;
							// CM Angles
							vcm = (m[11]/(m[11]+m[Atarget]))*sqrt(2*Ein/m[11]);
							vout[0] = sqrt(2*Efor[2][j]/m[Apart])*cos(theta);
							vout[1] = sqrt(2*Efor[2][j]/m[Apart])*sin(theta);
							theta_cm[2][j] = atan(vout[1]/(vout[0]-vcm));
							if (theta_cm[2][j] < 0) theta_cm[2][j] += pi;
      
							// Relativistic factors
							gamma[2][j] = (Ein + (m[11] + m[Atarget] - m[11+Atarget-Apart] - m[Apart]) - Eex[2][j] - Efor[2][j] + m[11+Atarget-Apart])/m[11+Atarget-Apart];
							beta[2][j] = sqrt(1. - 1./(gamma[2][j]*gamma[2][j]));
							
							
							if (strcmp(Histos[ih],"part")==0){ 
									CD_dEvsLabangle[(*FCD)[j].GetID()]->Fill(theta/pi*180.,(*FCD)[j].GetRingEnergy()[k]/1.e3);
									CD_dEvsLabangletotal->Fill(theta/pi*180.,(*FCD)[j].GetRingEnergy()[k]/1.e3);
							
							
									if (((*FCD)[j].GetID()==2)||((*FCD)[j].GetID()==3)){
										All_dEvsLabangletotal->Fill(theta,(*FCD)[j].GetRingEnergy()[k]);
//										NPT_Eextotal->Fill(Eex[2][j]/1.e3);
										
										for (Int_t iAngle=0; iAngle<20.;iAngle++)
										{
//											if ((theta_cm[2][j]>(iAngle*3./180.*pi + 85./180.*pi)) && (theta_cm[2][j]<((iAngle+1)*3./180.*pi + 85./180.*pi))){
//												  NPT_Eex_Angle[iAngle]->Fill(Eex[2][j]/1.e3);
//											}
											if ((theta>(iAngle*1./180.*pi + 8./180.*pi)) && (theta<((iAngle+1)*1./180.*pi + 8./180.*pi)) && (Apart==11)){
												Beryllium_Eex_Angle_Stop_CD[iAngle]->Fill(Eex[2][j]/1.e3);
											}
										}
										if (((*FCD)[j].GetRingEnergy()[k] > Be_cut[2]->Eval(theta)) && ((*FCD)[j].GetRingEnergy()[k] < Be_cut[3]->Eval(theta))) Be11_Angle->Fill(theta/pi*180.);
									}
									if (sim) {
										for (Int_t ie=0;ie<5;ie++){
											// Protons
											if ((par[0][ie]==1) || (par[0][ie]==3)){
												Q = m[11] + m[Atarget] - m[Apart] - m[11+Atarget-Apart];
												Pp = sqrt(2*m[Apart]*Efor[0][ie]);
												Pin = sqrt(2*m[11]*Ein);
												a = 1;
												b = 2*Pp*(cos(thetap[2][j])*cos(thetap[0][ie]) + sin(thetap[2][j])*sin(thetap[0][ie]));
												c = Pp*Pp - Pin*Pin;

												PN = (-b + sqrt(b*b - 4*a*c))/(2*a);
//												theta_coin = (sin(thetap[0][ie]) - cos(thetap[0][ie])*tan(thetap[2][j]))/(sin(thetain) - cos(thetain)*tan(thetap[2][j]));
//												PN = Pp*(theta_coin*cos(thetain) - cos(thetap[0][ie]))/cos(thetap[2][j]);
												EN = pow(PN,2)/(2*m[11+Atarget-Apart]);
//												Pin = Pp*theta_coin;
												Eex_coin = Ein + Q - Efor[0][ie] - EN;
//												Eex_coin = Ein_coin + Q - Efor[0][ie] - EN;
//											  	Eex_coin = Ein + Q - Efor[0][ie] - Efor[2][j];
												NPT_Eex_coin_Gamma_2700->Fill(Eex_coin/1.e3);
											}	
										}
									}
								}
								

								if (strcmp(Histos[ih],"coin")==0){	
									// Coincidence events
									for (Int_t ie=0;ie<5;ie++){
										// Protons
										if (par[0][ie]==1){
											Rafter = srim.RT[12]->Eval((*FCD)[j].GetStripEnergy()[l]);
											Efor[2][j] = srim.ET[12]->Eval(Rafter + ttarget/(2.*cos(theta)));
											Eex[2][j] = (pow(sqrt(2*m[12]*Efor[2][j])*cos(theta) + sqrt(2*m[1]*Efor[0][ie])*cos(thetap[0][ie]),2) + pow(sqrt(2*m[12]*Efor[2][j])*sin(theta) - sqrt(2*m[1]*Efor[0][ie])*sin(thetap[0][ie]),2))/(2*m[11]) - Efor[2][j] - Efor[0][ie] + (m[11] + m[2] - m[12] - m[1]);
											Proton_Eex_coin[ID[0][ie]]->Fill(Eex[2][j]/1.e3);
											Proton_Eex_cointotal->Fill(Eex[2][j]/1.e3);
											Proton_EexvsLabangle_coin[ID[0][ie]]->Fill(theta/pi*180.,Eex[2][j]/1.e3);
		
											thetain_coin = atan((sqrt(2*m[12]*Efor[2][j])*sin(theta)*cos(phi) + sqrt(2*m[1]*Efor[0][ie])*sin(thetap[0][ie])*cos(phip[0][ie]))/(sqrt(2*m[12]*Efor[2][j])*cos(theta) + sqrt(2*m[1]*Efor[0][ie])*cos(thetap[0][ie])));
											Proton_ThetaIn_coin[ID[0][ie]]->Fill(thetain_coin/pi*180.);
											Proton_ThetaInvsEex_coin[ID[0][ie]]->Fill(thetain_coin/pi*180.,Eex[2][j]/1.e3);
											Proton_PhivsPhi_coin[ID[0][ie]]->Fill(phip[0][ie]/pi*180,phi/pi*180.);
//											Hitpattern_coin[ID[0][ie]]->Fill(-63*tan(theta)*sin(phi),63*tan(theta)*cos(phi));
//											Hitpattern_cointotal->Fill(-63*tan(theta)*sin(phi),63*tan(theta)*cos(phi));
//											XY_coin[ID[0][ie]]->Fill(-(sin(phip[0][ie])+sin(phi))/2.,(cos(phip[0][ie])+cos(phi))/2.);
//											XY_cointotal->Fill(-(sin(phip[0][ie])+sin(phi))/2.,(cos(phip[0][ie])+cos(phi))/2.);
//											DetvsDet_coin->Fill(ID[0][ie],(*FCD)[j].GetID());	
										}
										// Deuterons
										if (par[0][ie]==2){
											Rafter = srim.RT[11]->Eval((*FCD)[j].GetStripEnergy()[l]);
											Efor[2][j] = srim.ET[11]->Eval(Rafter + ttarget/(2.*cos(theta)));
											Eex[2][j] = (pow(sqrt(2*m[11]*Efor[2][j])*cos(theta) + sqrt(2*m[2]*Efor[0][ie])*cos(thetap[0][ie]),2) + pow(sqrt(2*m[11]*Efor[2][j])*sin(theta) - sqrt(2*m[2]*Efor[0][ie])*sin(thetap[0][ie]),2))/(2*m[11]) - Efor[2][j] - Efor[0][ie] + (m[11] + m[2] - m[11] - m[2]);
											Deuteron_Eex_coin[ID[0][ie]]->Fill(Eex[2][j]/1.e3);
											Deuteron_Eex_cointotal->Fill(Eex[2][j]/1.e3);
											Deuteron_EexvsLabangle_coin[ID[0][ie]]->Fill(theta/pi*180.,Eex[2][j]/1.e3);
											thetain_coin = atan((sqrt(2*m[11]*Efor[2][j])*sin(theta)*cos(phi) + sqrt(2*m[2]*Efor[0][ie])*sin(thetap[0][ie])*cos(phip[0][ie]))/(sqrt(2*m[11]*Efor[2][j])*cos(theta) + sqrt(2*m[2]*Efor[0][ie])*cos(thetap[0][ie])));
											Deuteron_ThetaIn_coin[ID[0][ie]]->Fill(thetain_coin/pi*180.);
											Deuteron_ThetaInvsEex_coin[ID[0][ie]]->Fill(thetain_coin/pi*180.,Eex[2][j]/1.e3);
											Deuteron_PhivsPhi_coin[ID[0][ie]]->Fill(phip[0][ie]/pi*180,phi/pi*180.);
//											Hitpattern_coin[ID[0][ie]]->Fill(-63*tan(theta)*sin(phi),63*tan(theta)*cos(phi));
//											Hitpattern_cointotal->Fill(-63*tan(theta)*sin(phi),63*tan(theta)*cos(phi));
//											XY_coin[ID[0][ie]]->Fill(-(sin(phip[0][ie])+sin(phi))/2.,(cos(phip[0][ie])+cos(phi))/2.);
//											XY_cointotal->Fill(-(sin(phip[0][ie])+sin(phi))/2.,(cos(phip[0][ie])+cos(phi))/2.);
//											DetvsDet_coin->Fill(ID[0][ie],(*FCD)[j].GetID());	
										}
										// Tritons
										if (par[0][ie]==3){
											Rafter = srim.RDL[10]->Eval((*FCD)[j].GetStripEnergy()[l]);
											Efor[2][j] = srim.EDL[10]->Eval(Rafter + .7/cos(theta));
											Rafter = srim.RT[10]->Eval(Efor[2][j]);
											Efor[2][j] = srim.ET[10]->Eval(Rafter + ttarget/(2.*cos(theta)));
											Px = sqrt(2*m[10]*Efor[2][j])*sin(theta)*cos(phi)+sqrt(2*m[3]*Efor[0][ie])*sin(thetap[0][ie])*cos(phip[0][ie]);
											Py = sqrt(2*m[10]*Efor[2][j])*sin(theta)*sin(phi)+sqrt(2*m[3]*Efor[0][ie])*sin(thetap[0][ie])*sin(phip[0][ie]);
											Pz = sqrt(2*m[10]*Efor[2][j])*cos(theta)+sqrt(2*m[3]*Efor[0][ie])*cos(thetap[0][ie]);
											Eex[2][j] = (pow(Px,2)+pow(Py,2)+pow(Pz,2))/(2*m[11]) - Efor[2][j] - Efor[0][ie] + (m[11] + m[2] - m[10] - m[3]);
											Ein_coin = (pow(Px,2)+pow(Py,2)+pow(Pz,2))/(2*m[11]);
											
											Triton_Eex_coin[ID[0][ie]]->Fill(Eex[2][j]/1.e3);
											Triton_Eex_cointotal->Fill(Eex[2][j]/1.e3);
											Triton_Ein1_cointotal->Fill(Ein_coin/1.e3);
											
											Ein_coin = Efor[2][j] + Efor[0][ie] - (m[11]+m[2]-m[10]-m[3]);
											Triton_Ein2_cointotal->Fill(Ein_coin/1.e3);
											
											Triton_EexvsLabangle_coin[ID[0][ie]]->Fill(theta/pi*180.,Eex[2][j]/1.e3);
											thetain_coin = atan(sqrt(pow(Px,2)+pow(Py,2))/Pz);
//											thetain_coin = atan((sqrt(2*m[10]*Efor[2][j])*sin(theta)*cos(phi) + sqrt(2*m[3]*Efor[0][ie])*sin(thetap[0][ie])*sin(phip[0][ie]))/(sqrt(2*m[10]*Efor[2][j])*cos(theta) + sqrt(2*m[3]*Efor[0][ie])*cos(thetap[0][ie])));
											if (Px>=0&&Py>=0) phiin_coin = atan(Py/Px);
											if (Px>=0&&Py<0) phiin_coin = 2*pi + atan(Py/Px);
											if (Px<0&&Py>=0) phiin_coin = pi + atan(Py/Px);
											if (Px<0&&Py<0) phiin_coin = pi + atan(Py/Px);
											
											Triton_ThetaIn_coin[ID[0][ie]]->Fill(thetain_coin/pi*180.);
											Triton_PhiIn_coin[ID[0][ie]]->Fill(phiin_coin/pi*180.);
											Triton_ThetaInvsEex_coin[ID[0][ie]]->Fill(thetain_coin/pi*180.,Eex[2][j]/1.e3);
											Triton_PhivsPhi_coin[ID[0][ie]]->Fill(phip[0][ie]/pi*180,phi/pi*180.);
											Hitpattern_coin[ID[0][ie]]->Fill(63*tan(theta)*cos(phi),63*tan(theta)*sin(phi));
											Hitpattern_cointotal->Fill(63*tan(theta)*cos(phi),63*tan(theta)*sin(phi));
											XY_coin[ID[0][ie]]->Fill((cos(phip[0][ie])+cos(phi))/2.,(sin(phip[0][ie])+sin(phi))/2.);
											XY_cointotal->Fill((cos(phip[0][ie])+cos(phi))/2.,(sin(phip[0][ie])+sin(phi))/2.);
											XY_Beam_coin[ID[0][ie]]->Fill(thetain_coin/pi*180.*cos(phiin_coin),thetain_coin/pi*180.*sin(phiin_coin));
											XY_Beam_cointotal->Fill(thetain_coin/pi*180.*cos(phiin_coin),thetain_coin/pi*180.*sin(phiin_coin));
											DetvsDet_coin->Fill(ID[0][ie],(*FCD)[j].GetID());	
										}
										// No punch through
										if (par[0][ie]==4){
//											Rafter = srim.RDL[11+Atarget-Apart]->Eval((*FCD)[j].GetStripEnergy()[l]);
//											Efor[2][j] = srim.EDL[11+Atarget-Apart]->Eval(Rafter + .7/cos(theta));
//											Rafter = srim.RT[11+Atarget-Apart]->Eval(Efor[2][j]);
//											Efor[2][j] = srim.ET[11+Atarget-Apart]->Eval(Rafter + ttarget/(2.*cos(theta)));
											
											Eex[2][j] = -Efor[2][j] - Efor[0][ie] + (m[11]+m[Atarget]-m[11+Atarget-Apart]-m[Apart]) + Ein;
											NPT_Eex_cointotal->Fill(Eex[0][ie]/1e3);
											NPT_EexvsEex1_cointotal->Fill(Eex[2][j]/1e3,Eex[0][ie]/1e3);
											NPT_EtotvsLabangletotal_coin->Fill(thetap[0][ie]/pi*180.,Efor[0][ie]/1e3);
											if (Apart!=4){
											Px = Pxin - (sqrt(2*m[Atarget+11-Apart-1]*Efor[2][j])*sin(theta)*cos(phi)+sqrt(2*m[Apart]*Efor[0][ie])*sin(thetap[0][ie])*cos(phip[0][ie]));
											Py = Pyin - (sqrt(2*m[Atarget+11-Apart-1]*Efor[2][j])*sin(theta)*sin(phi)+sqrt(2*m[Apart]*Efor[0][ie])*sin(thetap[0][ie])*sin(phip[0][ie]));
											Pz = Pzin - (sqrt(2*m[Atarget+11-Apart-1]*Efor[2][j])*cos(theta)+sqrt(2*m[Apart]*Efor[0][ie])*cos(thetap[0][ie]));
											
											Eex[2][j] = Ein - (pow(Px,2)+pow(Py,2)+pow(Pz,2))/(2*m[0]) - Efor[2][j] - Efor[0][ie] + (m[11] + m[Atarget] - m[11+Atarget-1-Apart] - m[Apart] - m[0]);
											}
//											cout << Eex[2][j] << endl;
											NPT_EexvsEex2_cointotal->Fill(Eex[2][j]/1e3,Eex[0][ie]/1e3);
										}
									}
								}
							}
	
							// Particle identification
							if ((*FCD)[j].GetEdet() > CDThreshold[(*FCD)[j].GetID()]){
								Ecor = (*FCD)[j].GetEdet() + (1. - cos(theta))*(*FCD)[j].GetStripEnergy()[l];
								dEcor = (*FCD)[j].GetStripEnergy()[l]*cos(theta);
								
								if ((strcmp(Histos[ih],"part")==0)||(strcmp(Histos[ih],"sie")==0)){ 
									CD_dEvsEcortotal->Fill(Ecor/1.e3,dEcor/1.e3);
									CD_dEvsEcortot[(*FCD)[j].GetID()]->Fill(Ecor/1.e3,dEcor/1.e3);
								}

								if ((dEcor > CD_cut[0][(*FCD)[j].GetID()]->Eval(Ecor)) && (dEcor < CD_cut[1][(*FCD)[j].GetID()]->Eval(Ecor))){
									par[2][j] = 1;
									Np++;
									Np_CD++;
	
									Rafter = srim.RT[1]->Eval((*FCD)[j].GetEdet() + (*FCD)[j].GetStripEnergy()[l]);
									Efor[2][j] = srim.ET[1]->Eval(Rafter + ttarget/(2.*cos(theta)));
									Px = sqrt(2*m[1]*Efor[2][j])*sin(theta)*cos(phi);
									Py = sqrt(2*m[1]*Efor[2][j])*sin(theta)*sin(phi);
									Pz = sqrt(2*m[1]*Efor[2][j])*cos(theta);
									Q = m[11] + m[Atarget] - m[1] - m[11+Atarget-1];
									EN = (pow(Pxin-Px,2) + pow(Pyin-Py,2) + pow(Pzin-Pz,2))/(2*m[11+Atarget-1]);
						
									Eex[2][j] = Ein + Q - Efor[2][j] - EN;
	
									// CM Angles
									vcm = (m[11]/(m[11]+m[Atarget]))*sqrt(2*Ein/m[11]);
									vout[0] = sqrt(2*Efor[2][j]/m[1])*cos(theta);
									vout[1] = sqrt(2*Efor[2][j]/m[1])*sin(theta);
									theta_cm[2][j] = atan(vout[1]/(vout[0]-vcm));
									if (theta_cm[2][j] < 0) theta_cm[2][j] += pi;
									
									gamma[2][j] = (Ein + (m[11] + m[Atarget] - m[11+Atarget-1] - m[1]) - Eex[0][j] - Efor[2][j] + m[11+Atarget-1])/m[11+Atarget-1];
									//gamma[2][j] = (Ein + (m[22] + m[2] - m[23] - m[1]) - Eex[0][j] - Efor[2][j] + m[23])/m[23];
									beta[2][j] = sqrt(1. - 1./(gamma[2][j]*gamma[2][j]));
									
									if (strcmp(Histos[ih],"part")==0){
										// Filling histograms
										CD_dEvsEcorp->Fill(Ecor/1.e3,dEcor/1.e3);
										Proton_EtotvsLabangle[(*FCD)[j].GetID()]->Fill(theta/pi*180.,Efor[2][j]/1.e3);
										Proton_EtotvsLabangletotal->Fill(theta/pi*180.,Efor[2][j]/1.e3);
										Proton_EtotvsCmangletotal->Fill(theta_cm[2][j]/pi*180.,Efor[0][j]/1.e3);
										Proton_Eex_CD[(*FCD)[j].GetID()]->Fill(Eex[2][j]/1.e3);
										Proton_Eextotal->Fill(Eex[2][j]/1.e3);
										Proton_EexvsLabangle_CD[(*FCD)[j].GetID()]->Fill(theta/pi*180.,Eex[2][j]/1.e3);
										/*
										for (Int_t iAngle=0; iAngle<20.;iAngle++){
											if ((theta_cm[2][j]>(iAngle*2.5/180.*pi + 50./180.*pi)) && (theta_cm[2][j]<((iAngle+1)*2.5/180.*pi + 50./180.*pi))) Proton_Eex_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
										}
										*/
									}
								}
								if ((dEcor > CD_cut[2][(*FCD)[j].GetID()]->Eval(Ecor)) && (dEcor < CD_cut[3][(*FCD)[j].GetID()]->Eval(Ecor))){
									par[2][j] = 2;
									Nd++;
									Nd_CD++;
									
									Rafter = srim.RT[1]->Eval((*FCD)[j].GetEdet() + (*FCD)[j].GetStripEnergy()[l]);
									Efor[2][j] = srim.ET[1]->Eval(Rafter + ttarget/(2.*cos(theta)));
									Px = sqrt(2*m[2]*Efor[2][j])*sin(theta)*cos(phi);
									Py = sqrt(2*m[2]*Efor[2][j])*sin(theta)*sin(phi);
									Pz = sqrt(2*m[2]*Efor[2][j])*cos(theta);
									Q = m[11] + m[Atarget] - m[2] - m[11+Atarget-2];
									EN = (pow(Pxin-Px,2) + pow(Pyin-Py,2) + pow(Pzin-Pz,2))/(2*m[11+Atarget-2]);
						
									Eex[2][j] = Ein + Q - Efor[2][j] - EN;
	
									// CM Angles
									vcm = (m[11]/(m[11]+m[Atarget]))*sqrt(2*Ein/m[11]);
									vout[0] = sqrt(2*Efor[2][j]/m[2])*cos(theta);
									vout[1] = sqrt(2*Efor[2][j]/m[2])*sin(theta);
									theta_cm[2][j] = atan(vout[1]/(vout[0]-vcm));
									if (theta_cm[2][j] < 0) theta_cm[2][j] += pi;
									
									gamma[2][j] = (Ein + (m[11] + m[Atarget] - m[11+Atarget-2] - m[2]) - Eex[0][j] - Efor[2][j] + m[11+Atarget-2])/m[11+Atarget-2];
									//gamma[2][j] = (Ein + (m[22] + m[2] - m[22] - m[2]) - Eex[0][j] - Efor[2][j] + m[22])/m[22];
									beta[2][j] = sqrt(1. - 1./(gamma[2][j]*gamma[2][j]));
									
									if (strcmp(Histos[ih],"part")==0){
										// Filling histograms
										CD_dEvsEcord->Fill(Ecor/1.e3,dEcor/1.e3);
										Deuteron_EtotvsLabangle[(*FCD)[j].GetID()]->Fill(theta/pi*180.,Efor[2][j]/1.e3);
										Deuteron_EtotvsLabangletotal->Fill(theta/pi*180.,Efor[2][j]/1.e3);
										Deuteron_EtotvsCmangletotal->Fill(theta_cm[2][j]/pi*180.,Efor[0][j]/1.e3);
										Deuteron_Eex_CD[(*FCD)[j].GetID()]->Fill(Eex[2][j]/1.e3);
										Deuteron_Eextotal->Fill(Eex[2][j]/1.e3);
										Deuteron_EexvsLabangle_CD[(*FCD)[j].GetID()]->Fill(theta/pi*180.,Eex[2][j]/1.e3);
										/*
										for (Int_t iAngle=0; iAngle<20.;iAngle++){
											if ((theta_cm[2][j]>(iAngle*2.5/180.*pi + 50./180.*pi)) && (theta_cm[2][j]<((iAngle+1)*2.5/180.*pi + 50./180.*pi))) Deuteron_Eex_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
										}
										*/
									}
								}
								if ((dEcor > CD_cut[4][(*FCD)[j].GetID()]->Eval(Ecor)) && (dEcor < CD_cut[5][(*FCD)[j].GetID()]->Eval(Ecor))){
									par[2][j] = 3;
									Nt++;
									Nt_CD++;
									
									Rafter = srim.RT[1]->Eval((*FCD)[j].GetEdet() + (*FCD)[j].GetStripEnergy()[l]);
									Efor[2][j] = srim.ET[1]->Eval(Rafter + ttarget/(2.*cos(theta)));
									Px = sqrt(2*m[3]*Efor[2][j])*sin(theta)*cos(phi);
									Py = sqrt(2*m[3]*Efor[2][j])*sin(theta)*sin(phi);
									Pz = sqrt(2*m[3]*Efor[2][j])*cos(theta);
									Q = m[11] + m[Atarget] - m[3] - m[11+Atarget-3];
									EN = (pow(Pxin-Px,2) + pow(Pyin-Py,2) + pow(Pzin-Pz,2))/(2*m[11+Atarget-3]);
						
									Eex[2][j] = Ein + Q - Efor[2][j] - EN;
	
									// CM Angles
									vcm = (m[11]/(m[11]+m[Atarget]))*sqrt(2*Ein/m[11]);
									vout[0] = sqrt(2*Efor[2][j]/m[3])*cos(theta);
									vout[1] = sqrt(2*Efor[2][j]/m[3])*sin(theta);
									theta_cm[2][j] = atan(vout[1]/(vout[0]-vcm));
									if (theta_cm[2][j] < 0) theta_cm[2][j] += pi;
									
									gamma[2][j] = (Ein + (m[11] + m[Atarget] - m[11+Atarget-3] - m[3]) - Eex[0][j] - Efor[2][j] + m[11+Atarget-3])/m[11+Atarget-3];
									beta[2][j] = sqrt(1. - 1./(gamma[2][j]*gamma[2][j]));
									
									if (strcmp(Histos[ih],"part")==0){
										// Filling histograms
										CD_dEvsEcort->Fill(Ecor/1.e3,dEcor/1.e3);
										Triton_EtotvsLabangle[(*FCD)[j].GetID()]->Fill(theta/pi*180.,Efor[2][j]/1.e3);
										Triton_EtotvsLabangletotal->Fill(theta/pi*180.,Efor[2][j]/1.e3);
										Triton_EtotvsCmangletotal->Fill(theta_cm[2][j]/pi*180.,Efor[0][j]/1.e3);
										Triton_Eex_CD[(*FCD)[j].GetID()]->Fill(Eex[2][j]/1.e3);
										Triton_Eextotal->Fill(Eex[2][j]/1.e3);
										Triton_EexvsLabangle_CD[(*FCD)[j].GetID()]->Fill(theta/pi*180.,Eex[2][j]/1.e3);
										/*
										for (Int_t iAngle=0; iAngle<20.;iAngle++){
											if ((theta_cm[2][j]>(iAngle*2.5/180.*pi + 50./180.*pi)) && (theta_cm[2][j]<((iAngle+1)*2.5/180.*pi + 50./180.*pi))) Triton_Eex_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
										}
										*/
									}
								}
							}
						}
					}
				}
			}
		}
		// End cd


		// Gammas
		for(unsigned int j=0; j<Miniball->size(); j++){
			if((*Miniball)[j].GetCluID()<0 || (*Miniball)[j].GetCluID()>7){
				cout << "Error in entry "<<i<<": "<<j<<". Germanium detector Cluster id is wrong: "<<(*Miniball)[j].GetCluID()<<endl;
				continue;
			}

			// Filling histograms
			for(unsigned int k=0; k<(*Miniball)[j].GetCrystal().size();k++){
				if( (*Miniball)[j].GetCrystal()[k].GetCryID() <0 || (*Miniball)[j].GetCrystal()[k].GetCryID() >2){
					cout << "Error in entry "<<i<<": "<<j<<". Germanium detector Crystal id is wrong: "<<(*Miniball)[j].GetCrystal()[k].GetCryID()<<endl;
					continue;
				}
				
				for (Int_t ih=0;ih<Histos.size();ih++){
					if (strcmp(Histos[ih],"gee")==0){
						Gamma_Etotal->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
						Gamma_ECrystals1[(*Miniball)[j].GetCluID()][(*Miniball)[j].GetCrystal()[k].GetCryID()]->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
						Gamma_ECrystals2[(*Miniball)[j].GetCluID()][(*Miniball)[j].GetCrystal()[k].GetCryID()]->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
					}
	
					for (Int_t l=0;l<1;l++){
						// Sometimes segment ID is not set. Why? Who knows. It doesn't
						// matter much for the analysis, since it only influences the lab
						// angle, and at most 1--2 degrees with current (IS430) setup.
						// Events could have been skipped, but for the current analysis it
						// is better to have a small uncertainty in the angle than drop the
						// event. If dropping is wanted, use the for condition:
						//
						//   for (Int_t l = 0; l < (*Miniball)[j].GetCrystal()[k].GetSegID().size(); l++) {
						//
						// above.
						//
						// If the segment ID is not set from the data file, then it is
						// assumed to be a truly random number within 0--5 here
						// (<http://xkcd.com/221/>).
						int pseudoSegID;
						if (0 == (*Miniball)[j].GetCrystal()[k].GetSegID().size()) {
							pseudoSegID = 4;
						} else {
							pseudoSegID = (*Miniball)[j].GetCrystal()[k].GetSegID()[l];
						}

						if (strcmp(Histos[ih],"gee")==0) Gamma_ESegments[(*Miniball)[j].GetCluID()][(*Miniball)[j].GetCrystal()[k].GetCryID()][pseudoSegID]->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
	
						// Gammas in coincidence with particles
						for (Int_t idet=0;idet<3;idet++){
							for (Int_t ie=0;ie<5;ie++){
								if ((par[idet][ie]==1)&&(ID[idet][ie]!=0)){
									if (strcmp(Histos[ih],"part")==0){
										EGamDS = gamma[idet][ie]*(*Miniball)[j].GetCrystal()[k].GetCore()*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[j].GetCluID()][(*Miniball)[j].GetCrystal()[k].GetCryID()][pseudoSegID]));
										Proton_ESegments[(*Miniball)[j].GetCluID()][(*Miniball)[j].GetCrystal()[k].GetCryID()][pseudoSegID]->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
										Proton_EGamma->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
										Proton_EDS->Fill(EGamDS/1.e3);
										Gamma_EDS->Fill(EGamDS/1.e3);
										Proton_EexvsEGamma->Fill(Eex[idet][ie]/1.e3,(*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
										Proton_EexvsEDS->Fill(Eex[idet][ie]/1.e3,EGamDS/1.e3);
										// Gamma Gamma coincidence
										if ((Miniball->size()==2) && ((*Miniball)[j].GetCrystal().size()==1) && (j==0) && (k==0) && (l==0)){
											int pseudoSegID2;
											if (0 == (*Miniball)[1].GetCrystal()[0].GetSegID().size()) {
												pseudoSegID2 = 4;
											} else {
												pseudoSegID2 = (*Miniball)[1].GetCrystal()[0].GetSegID()[0];
											}

											EGamDS_2 = gamma[idet][ie]*(*Miniball)[1].GetCrystal()[0].GetCore()*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[1].GetCluID()][(*Miniball)[1].GetCrystal()[0].GetCryID()][pseudoSegID2]));
											EGamma_total = (*Miniball)[0].GetCrystal()[0].GetCore() + (*Miniball)[1].GetCrystal()[0].GetCore();
											EGamDS_total = gamma[idet][ie]*EGamma_total*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[0].GetCluID()][(*Miniball)[0].GetCrystal()[0].GetCryID()][(*Miniball)[0].GetCrystal()[0].GetSegID()[0]]));
											Proton_EGammavsEGamma->Fill((*Miniball)[0].GetCrystal()[0].GetCore()/1.e3,(*Miniball)[1].GetCrystal()[0].GetCore()/1.e3);
											Proton_EDSvsEDS->Fill(EGamDS/1.e3,EGamDS_2/1.e3);
											Proton_EGamma_Total->Fill(EGamma_total/1.e3);
											Proton_EDS_Total->Fill(EGamDS_total/1.e3);
										}
										else if ((Miniball->size()==1) && ((*Miniball)[j].GetCrystal().size()==2) && (j==0) && (k==0) & (l==0)){
											EGamDS_2 = gamma[idet][ie]*(*Miniball)[0].GetCrystal()[1].GetCore()*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[0].GetCluID()][(*Miniball)[0].GetCrystal()[1].GetCryID()][(*Miniball)[0].GetCrystal()[1].GetSegID()[0]]));
											EGamma_total = (*Miniball)[0].GetCrystal()[0].GetCore() + (*Miniball)[0].GetCrystal()[1].GetCore();
											EGamDS_total = gamma[idet][ie]*EGamma_total*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[0].GetCluID()][(*Miniball)[0].GetCrystal()[0].GetCryID()][(*Miniball)[0].GetCrystal()[0].GetSegID()[0]]));
											Proton_EGammavsEGamma->Fill((*Miniball)[0].GetCrystal()[0].GetCore()/1.e3,(*Miniball)[0].GetCrystal()[1].GetCore()/1.e3);
											Proton_EDSvsEDS->Fill(EGamDS/1.e3,EGamDS_2/1.e3);
											Proton_EGamma_Total->Fill(EGamma_total/1.e3);
											Proton_EDS_Total->Fill(EGamDS_total/1.e3);
										}
										if ((Miniball->size()==1) && ((*Miniball)[j].GetCrystal().size()==1)){
											Proton_EGamma_Total->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
											Proton_EGamma_1->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
											Proton_EDS_Total->Fill(EGamDS/1.e3);
									 		Proton_EDS_1->Fill(EGamDS/1.e3);
											Proton_EexvsEGamma_1->Fill(Eex[idet][ie]/1.e3,(*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
											Proton_EexvsEDS_1->Fill(Eex[idet][ie]/1.e3,EGamDS/1.e3);
											if ((EGamDS<2160)&&(2060<EGamDS)){
												Proton_Eex_Gamma_2100_1->Fill(Eex[idet][ie]/1.e3);
											}
											if ((EGamDS<2790)&&(2690<EGamDS)){
												Proton_Eex_Gamma_2700_1->Fill(Eex[idet][ie]/1.e3);
											}
										}
										else if (((Miniball->size()==2) && ((*Miniball)[0].GetCrystal().size()==1) && ((*Miniball)[1].GetCrystal().size()==1)) || ((Miniball->size()==1) && (((*Miniball)[0].GetCrystal().size()==2)))){
											Proton_EGamma_2->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
									 		Proton_EDS_2->Fill(EGamDS/1.e3);
										}
										else{
											Proton_EGamma_3->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
									 		Proton_EDS_3->Fill(EGamDS/1.e3);
										}

									
									if ((MBTheta[(*Miniball)[0].GetCluID()][(*Miniball)[0].GetCrystal()[0].GetCryID()][(*Miniball)[0].GetCrystal()[0].GetSegID()[0]] > -10./180.*pi) && (MBTheta[(*Miniball)[0].GetCluID()][(*Miniball)[0].GetCrystal()[0].GetCryID()][(*Miniball)[0].GetCrystal()[0].GetSegID()[0]] < 360./180.*pi)){
										if (((*Miniball)[j].GetCrystal()[k].GetCore()<520)&&(490<(*Miniball)[j].GetCrystal()[k].GetCore())){
											Proton_Eex_Gamma_511->Fill(Eex[idet][ie]/1.e3);
											for (Int_t iAngle=0; iAngle<20.;iAngle++)
											{
												if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 50./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 50./180.*pi))) Proton_Eex_511_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
											}
										}
										if ((((*Miniball)[j].GetCrystal()[k].GetCore()<535)&&(520<(*Miniball)[j].GetCrystal()[k].GetCore()))||(((*Miniball)[j].GetCrystal()[k].GetCore()<490)&&(475<(*Miniball)[j].GetCrystal()[k].GetCore()))){
											Proton_Eex_Gamma_511_Background->Fill(Eex[idet][ie]/1.e3);
											for (Int_t iAngle=0; iAngle<20.;iAngle++)
											{
												if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 50./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 50./180.*pi))) Proton_Eex_511_Angle_Background[iAngle]->Fill(Eex[idet][ie]/1.e3);
											}
										}
										if ((EGamDS<2160)&&(2060<EGamDS)){
											Proton_Eex_Gamma_2100->Fill(Eex[idet][ie]/1.e3);
											Proton_EtotvsCmangle_Gamma_2100->Fill(theta_cm[idet][ie]/pi*180.,Efor[idet][ie]/1.e3);
											Proton_MBTheta_2100->Fill(MBTheta[(*Miniball)[0].GetCluID()][(*Miniball)[0].GetCrystal()[0].GetCryID()][(*Miniball)[0].GetCrystal()[0].GetSegID()[0]]/pi*180.);

											for (Int_t iAngle=0; iAngle<20.;iAngle++)
											{
												if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 50./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 50./180.*pi))) Proton_Eex_2100_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
											}
										}
										if ((EGamDS<2260)&&(2160<EGamDS)){
											for (Int_t iAngle=0; iAngle<20.;iAngle++)
											{
												if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 50./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 50./180.*pi))) Proton_Eex_2100_Angle_Background[iAngle]->Fill(Eex[idet][ie]/1.e3);
											}
										}
										if ((EGamDS<2790)&&(2690<EGamDS)){
											Proton_Eex_Gamma_2700->Fill(Eex[idet][ie]/1.e3);
											Proton_EtotvsCmangle_Gamma_2700->Fill(theta_cm[idet][ie]/pi*180.,Efor[idet][ie]/1.e3);
											Proton_MBTheta_2700->Fill(MBTheta[(*Miniball)[0].GetCluID()][(*Miniball)[0].GetCrystal()[0].GetCryID()][(*Miniball)[0].GetCrystal()[0].GetSegID()[0]]/pi*180.);
											for (Int_t iAngle=0; iAngle<20.;iAngle++)
											{
												if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 50./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 50./180.*pi))) Proton_Eex_2700_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
											}
										}
									}
									
									}
									
									if (strcmp(Histos[ih],"time")==0){
										Gamma_EvsdtProton->Fill(25*((*Miniball)[j].GetCrystal()[k].GetTime() - time[idet][ie]),(*Miniball)[j].GetCrystal()[k].GetCore()/1000.);
										if ((((*Miniball)[j].GetCrystal()[k].GetCore()<520)&&(500<(*Miniball)[j].GetCrystal()[k].GetCore()))&&(1800<Eex[idet][ie])) Proton_dt_511->Fill(25*((*Miniball)[j].GetCrystal()[k].GetTime() - time[idet][ie]));
									}
								}
								else if ((par[idet][ie]==2)&&(ID[idet][ie]!=0)){
									if (strcmp(Histos[ih],"part")==0){
										EGamDS = gamma[idet][ie]*(*Miniball)[j].GetCrystal()[k].GetCore()*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[j].GetCluID()][(*Miniball)[j].GetCrystal()[k].GetCryID()][pseudoSegID]));
										Deuteron_ESegments[(*Miniball)[j].GetCluID()][(*Miniball)[j].GetCrystal()[k].GetCryID()][pseudoSegID]->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
										Deuteron_EGamma->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
										Deuteron_EDS->Fill(EGamDS/1.e3);
										Gamma_EDS->Fill(EGamDS/1.e3);
										Deuteron_EexvsEGamma->Fill(Eex[idet][ie]/1.e3,(*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
										Deuteron_EexvsEDS->Fill(Eex[idet][ie]/1.e3,EGamDS/1.e3);
										// Gamma Gamma coincidence
										if ((Miniball->size()==2) && ((*Miniball)[j].GetCrystal().size()==1) && (j==0) && (k==0) && (l==0)){
											EGamDS_2 = gamma[idet][ie]*(*Miniball)[1].GetCrystal()[0].GetCore()*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[1].GetCluID()][(*Miniball)[1].GetCrystal()[0].GetCryID()][(*Miniball)[1].GetCrystal()[0].GetSegID()[0]]));
											EGamma_total = (*Miniball)[0].GetCrystal()[0].GetCore() + (*Miniball)[1].GetCrystal()[0].GetCore();
											EGamDS_total = gamma[idet][ie]*EGamma_total*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[0].GetCluID()][(*Miniball)[0].GetCrystal()[0].GetCryID()][(*Miniball)[0].GetCrystal()[0].GetSegID()[0]]));
											Deuteron_EGammavsEGamma->Fill((*Miniball)[0].GetCrystal()[0].GetCore()/1.e3,(*Miniball)[1].GetCrystal()[0].GetCore()/1.e3);
											Deuteron_EDSvsEDS->Fill(EGamDS/1.e3,EGamDS_2/1.e3);
											Deuteron_EGamma_Total->Fill(EGamma_total/1.e3);
											Deuteron_EDS_Total->Fill(EGamDS_total/1.e3);
										}
										else if ((Miniball->size()==1) && ((*Miniball)[j].GetCrystal().size()==2) && (j==0) && (k==0) & (l==0)){
											EGamDS_2 = gamma[idet][ie]*(*Miniball)[0].GetCrystal()[1].GetCore()*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[0].GetCluID()][(*Miniball)[0].GetCrystal()[1].GetCryID()][(*Miniball)[0].GetCrystal()[1].GetSegID()[0]]));
											EGamma_total = (*Miniball)[0].GetCrystal()[0].GetCore() + (*Miniball)[0].GetCrystal()[1].GetCore();
											EGamDS_total = gamma[idet][ie]*EGamma_total*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[0].GetCluID()][(*Miniball)[0].GetCrystal()[0].GetCryID()][(*Miniball)[0].GetCrystal()[0].GetSegID()[0]]));
											Deuteron_EGammavsEGamma->Fill((*Miniball)[0].GetCrystal()[0].GetCore()/1.e3,(*Miniball)[0].GetCrystal()[1].GetCore()/1.e3);
											Deuteron_EDSvsEDS->Fill(EGamDS/1.e3,EGamDS_2/1.e3);
											Deuteron_EGamma_Total->Fill(EGamma_total/1.e3);
											Deuteron_EDS_Total->Fill(EGamDS_total/1.e3);
										}
										if ((Miniball->size()==1) && ((*Miniball)[j].GetCrystal().size()==1)){
											Deuteron_EGamma_Total->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
											Deuteron_EGamma_1->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
											Deuteron_EDS_Total->Fill(EGamDS/1.e3);
									 		Deuteron_EDS_1->Fill(EGamDS/1.e3);
											Deuteron_EexvsEGamma_1->Fill(Eex[idet][ie]/1.e3,(*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
											Deuteron_EexvsEDS_1->Fill(Eex[idet][ie]/1.e3,EGamDS/1.e3);
											if ((EGamDS<340)&&(300<EGamDS)){
												Deuteron_Eex_Gamma_320_1->Fill(Eex[idet][ie]/1.e3);
											}
										}
										else if (((Miniball->size()==2) && ((*Miniball)[0].GetCrystal().size()==1) && ((*Miniball)[1].GetCrystal().size()==1)) || ((Miniball->size()==1) && (((*Miniball)[0].GetCrystal().size()==2)))){
											Deuteron_EGamma_2->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
									 		Deuteron_EDS_2->Fill(EGamDS/1.e3);
										}
										else{
											Deuteron_EGamma_3->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
									 		Deuteron_EDS_3->Fill(EGamDS/1.e3);
										}

										
									if ((MBTheta[(*Miniball)[0].GetCluID()][(*Miniball)[0].GetCrystal()[0].GetCryID()][(*Miniball)[0].GetCrystal()[0].GetSegID()[0]] > -10./180.*pi) && (MBTheta[(*Miniball)[0].GetCluID()][(*Miniball)[0].GetCrystal()[0].GetCryID()][(*Miniball)[0].GetCrystal()[0].GetSegID()[0]] < 360./180.*pi)){
										if ((EGamDS<340)&&(300<EGamDS)){
											Deuteron_Eex_Gamma_320->Fill(Eex[idet][ie]/1.e3);
											Deuteron_EtotvsCmangle_Gamma_320->Fill(theta_cm[idet][ie]/pi*180.,Efor[idet][ie]/1.e3);
											Deuteron_MBTheta_320->Fill(MBTheta[(*Miniball)[0].GetCluID()][(*Miniball)[0].GetCrystal()[0].GetCryID()][(*Miniball)[0].GetCrystal()[0].GetSegID()[0]]/pi*180.);
											for (Int_t iAngle=0; iAngle<20.;iAngle++)
											{
												if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 60./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 60./180.*pi))) Deuteron_Eex_320_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
											}
										}
										if ((EGamDS<380)&&(340<EGamDS)){
											for (Int_t iAngle=0; iAngle<20.;iAngle++)
											{
												if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 60./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 60./180.*pi))) Deuteron_Eex_320_Angle_Background[iAngle]->Fill(Eex[idet][ie]/1.e3);
											}
										}
									}
									
									}
									if (strcmp(Histos[ih],"time")==0){
										Gamma_EvsdtDeuteron->Fill(25*((*Miniball)[j].GetCrystal()[k].GetTime() - time[idet][ie]),(*Miniball)[j].GetCrystal()[k].GetCore()/1000.);
									}
								}
								else if ((par[idet][ie]==3)&&(ID[idet][ie]!=0)){
									if (strcmp(Histos[ih],"part")==0){
										EGamDS = gamma[idet][ie]*(*Miniball)[j].GetCrystal()[k].GetCore()*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[j].GetCluID()][(*Miniball)[j].GetCrystal()[k].GetCryID()][pseudoSegID]));
										Triton_EGamma->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
										Triton_EDS->Fill(EGamDS/1.e3);
										Gamma_EDS->Fill(EGamDS/1.e3);
										Triton_EexvsEGamma->Fill(Eex[idet][ie]/1.e3,(*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
										Triton_EexvsEDS->Fill(Eex[idet][ie]/1.e3,EGamDS/1.e3);
										// Gamma Gamma coincidence
										if ((Miniball->size()==2) && ((*Miniball)[j].GetCrystal().size()==1) && (j==0) && (k==0) && (l==0)){
											EGamDS_2 = gamma[idet][ie]*(*Miniball)[1].GetCrystal()[0].GetCore()*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[1].GetCluID()][(*Miniball)[1].GetCrystal()[0].GetCryID()][(*Miniball)[1].GetCrystal()[0].GetSegID()[0]]));
											EGamma_total = (*Miniball)[0].GetCrystal()[0].GetCore() + (*Miniball)[1].GetCrystal()[0].GetCore();
											EGamDS_total = gamma[idet][ie]*EGamma_total*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[0].GetCluID()][(*Miniball)[0].GetCrystal()[0].GetCryID()][(*Miniball)[0].GetCrystal()[0].GetSegID()[0]]));
											Triton_EGammavsEGamma->Fill((*Miniball)[0].GetCrystal()[0].GetCore()/1.e3,(*Miniball)[1].GetCrystal()[0].GetCore()/1.e3);
											Triton_EDSvsEDS->Fill(EGamDS/1.e3,EGamDS_2/1.e3);
											Triton_EGamma_Total->Fill(EGamma_total/1.e3);
											Triton_EDS_Total->Fill(EGamDS_total/1.e3);
											
											if ((3300<EGamDS) && (EGamDS<3500)) Triton_EDS_coin_3300->Fill(EGamDS_2/1.e3);
											if ((3300<EGamDS_2) && (EGamDS_2<3500)) Triton_EDS_coin_3300->Fill(EGamDS/1.e3);
										}
										else if ((Miniball->size()==1) && ((*Miniball)[j].GetCrystal().size()==2) && (j==0) && (k==0) & (l==0)){
											EGamDS_2 = gamma[idet][ie]*(*Miniball)[0].GetCrystal()[1].GetCore()*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[0].GetCluID()][(*Miniball)[0].GetCrystal()[1].GetCryID()][(*Miniball)[0].GetCrystal()[1].GetSegID()[0]]));
											EGamma_total = (*Miniball)[0].GetCrystal()[0].GetCore() + (*Miniball)[0].GetCrystal()[1].GetCore();
											EGamDS_total = gamma[idet][ie]*EGamma_total*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[0].GetCluID()][(*Miniball)[0].GetCrystal()[0].GetCryID()][(*Miniball)[0].GetCrystal()[0].GetSegID()[0]]));
											Triton_EGammavsEGamma->Fill((*Miniball)[0].GetCrystal()[0].GetCore()/1.e3,(*Miniball)[0].GetCrystal()[1].GetCore()/1.e3);
											Triton_EDSvsEDS->Fill(EGamDS/1.e3,EGamDS_2/1.e3);
											Triton_EGamma_Total->Fill(EGamma_total/1.e3);
											Triton_EDS_Total->Fill(EGamDS_total/1.e3);
											
											if ((3300<EGamDS) && (EGamDS<3500)) Triton_EDS_coin_3300->Fill(EGamDS_2/1.e3);
											if ((3300<EGamDS_2) && (EGamDS_2<3500)) Triton_EDS_coin_3300->Fill(EGamDS/1.e3);
										}
										if ((Miniball->size()==1) && ((*Miniball)[j].GetCrystal().size()==1)){
											Triton_EGamma_Total->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
											Triton_EGamma_1->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
											Triton_EDS_Total->Fill(EGamDS/1.e3);
									 		Triton_EDS_1->Fill(EGamDS/1.e3);
											Triton_EexvsEGamma_1->Fill(Eex[idet][ie]/1.e3,(*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
											Triton_EexvsEDS_1->Fill(Eex[idet][ie]/1.e3,EGamDS/1.e3);
											if ((EGamDS<3500)&&(3300<EGamDS)){
												Triton_Eex_Gamma_3300_1->Fill(Eex[idet][ie]/1.e3);
											}
											if ((EGamDS<6200)&&(4400<EGamDS)){
												Triton_Eex_Gamma_6000_comp_1->Fill(Eex[idet][ie]/1.e3);
											}
										}
										else if (((Miniball->size()==2) && ((*Miniball)[0].GetCrystal().size()==1) && ((*Miniball)[1].GetCrystal().size()==1)) || ((Miniball->size()==1) && (((*Miniball)[0].GetCrystal().size()==2)))){
											Triton_EGamma_2->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
									 		Triton_EDS_2->Fill(EGamDS/1.e3);
										}
										else{
											Triton_EGamma_3->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
									 		Triton_EDS_3->Fill(EGamDS/1.e3);
										}
									
									if ((MBTheta[(*Miniball)[0].GetCluID()][(*Miniball)[0].GetCrystal()[0].GetCryID()][(*Miniball)[0].GetCrystal()[0].GetSegID()[0]] > -10./180.*pi) && (MBTheta[(*Miniball)[0].GetCluID()][(*Miniball)[0].GetCrystal()[0].GetCryID()][(*Miniball)[0].GetCrystal()[0].GetSegID()[0]] < 360./180.*pi)){
										if ((EGamDS<3500)&&(3300<EGamDS)){
											Triton_MBTheta_3300->Fill(MBTheta[(*Miniball)[0].GetCluID()][(*Miniball)[0].GetCrystal()[0].GetCryID()][(*Miniball)[0].GetCrystal()[0].GetSegID()[0]]/pi*180.);
												Triton_Eex_Gamma_3300->Fill(Eex[idet][ie]/1.e3);
												Triton_EtotvsCmangle_Gamma_3300->Fill(theta_cm[idet][ie]/pi*180.,Efor[idet][ie]/1.e3);
											for (Int_t iAngle=0; iAngle<20.;iAngle++)
											{
												if ((theta_cm[0][j]>(iAngle*3./180.*pi + 60./180.*pi)) && (theta_cm[0][j]<((iAngle+1)*3./180.*pi + 60./180.*pi))) Triton_Eex_3300_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
											}
										}
										if ((EGamDS<3500)&&(2200<EGamDS)){
											Triton_Eex_Gamma_3300_comp->Fill(Eex[idet][ie]/1.e3);
											for (Int_t iAngle=0; iAngle<20.;iAngle++)
											{
												if ((theta_cm[0][j]>(iAngle*3./180.*pi + 60./180.*pi)) && (theta_cm[0][j]<((iAngle+1)*3./180.*pi + 60./180.*pi))) Triton_Eex_3300_comp_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
											}
										}
										if ((EGamDS<3000)&&(2800<EGamDS)){
											for (Int_t iAngle=0; iAngle<20.;iAngle++)
											{
												if ((theta_cm[0][j]>(iAngle*3./180.*pi + 60./180.*pi)) && (theta_cm[0][j]<((iAngle+1)*3./180.*pi + 60./180.*pi))) Triton_Eex_2900_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
											}
										}
										if ((EGamDS<2750)&&(2500<EGamDS)){
											for (Int_t iAngle=0; iAngle<20.;iAngle++)
											{
												if ((theta_cm[0][j]>(iAngle*3./180.*pi + 60./180.*pi)) && (theta_cm[0][j]<((iAngle+1)*3./180.*pi + 60./180.*pi))) Triton_Eex_2600_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
											}
										}
										if ((EGamDS<6200)&&(5800<EGamDS)){
											for (Int_t iAngle=0; iAngle<20.;iAngle++)
											{
												if ((theta_cm[0][j]>(iAngle*3./180.*pi + 60./180.*pi)) && (theta_cm[0][j]<((iAngle+1)*3./180.*pi + 60./180.*pi))) Triton_Eex_6000_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
											}
										}
										if ((EGamDS<6200)&&(4400<EGamDS)){
											Triton_Eex_Gamma_6000_comp->Fill(Eex[idet][ie]/1.e3);
											Triton_EtotvsCmangle_Gamma_6000->Fill(theta_cm[idet][ie]/pi*180.,Efor[idet][ie]/1.e3);
											for (Int_t iAngle=0; iAngle<20.;iAngle++)
											{
												if ((theta_cm[0][j]>(iAngle*3./180.*pi + 60./180.*pi)) && (theta_cm[0][j]<((iAngle+1)*3./180.*pi + 60./180.*pi))) Triton_Eex_6000_comp_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
											}
										}
									}

									}
									if (strcmp(Histos[ih],"time")==0){
										Gamma_EvsdtTriton->Fill(25*((*Miniball)[j].GetCrystal()[k].GetTime() - time[idet][ie]),(*Miniball)[j].GetCrystal()[k].GetCore()/1000.);
									}
								}
								else if ((par[idet][ie]==4)&&(ID[idet][ie]!=0)){
									EGamDS = gamma[idet][ie]*(*Miniball)[j].GetCrystal()[k].GetCore()*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[j].GetCluID()][(*Miniball)[j].GetCrystal()[k].GetCryID()][pseudoSegID]));
									if (strcmp(Histos[ih],"part")==0){
										if ((idet==1) && ((ID[idet][ie]==2) || (ID[idet][ie]==3))){
											NPT_EGamma_Back->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
											NPT_EDS_Back->Fill(EGamDS/1.e3);
											NPT_EexvsEGamma_Back->Fill(Eex[idet][ie]/1.e3,(*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
											NPT_EexvsEDS_Back->Fill(Eex[idet][ie]/1.e3,EGamDS/1.e3);
										}
										else if (idet==0){
											NPT_EGamma->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
											NPT_EDS->Fill(EGamDS/1.e3);
											NPT_EGamma_Low->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
											NPT_EDS_Low->Fill(EGamDS/1.e3);
											NPT_EGamma_High->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
											NPT_EDS_High->Fill(EGamDS/1.e3);
											NPT_EexvsEGamma->Fill(Eex[idet][ie]/1.e3,(*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
											NPT_EexvsEDS->Fill(Eex[idet][ie]/1.e3,EGamDS/1.e3);
											// Gamma Gamma coincidence
											if ((Miniball->size()==2) && ((*Miniball)[j].GetCrystal().size()==1) && (j==0) && (k==0) && (l==0)){
												EGamDS_2 = gamma[idet][ie]*(*Miniball)[1].GetCrystal()[0].GetCore()*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[1].GetCluID()][(*Miniball)[1].GetCrystal()[0].GetCryID()][(*Miniball)[1].GetCrystal()[0].GetSegID()[0]]));
												EGamma_total = (*Miniball)[0].GetCrystal()[0].GetCore() + (*Miniball)[1].GetCrystal()[0].GetCore();
												EGamDS_total = gamma[idet][ie]*EGamma_total*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[0].GetCluID()][(*Miniball)[0].GetCrystal()[0].GetCryID()][(*Miniball)[0].GetCrystal()[0].GetSegID()[0]]));
												NPT_EGammavsEGamma->Fill((*Miniball)[0].GetCrystal()[0].GetCore()/1.e3,(*Miniball)[1].GetCrystal()[0].GetCore()/1.e3);
												NPT_EDSvsEDS->Fill(EGamDS/1.e3,EGamDS_2/1.e3);
												NPT_EDS_Total->Fill(EGamDS_total/1.e3);
											}
											else if ((Miniball->size()==1) && ((*Miniball)[j].GetCrystal().size()==2) && (j==0) && (k==0) & (l==0)){
												EGamDS_2 = gamma[idet][ie]*(*Miniball)[0].GetCrystal()[1].GetCore()*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[0].GetCluID()][(*Miniball)[0].GetCrystal()[1].GetCryID()][(*Miniball)[0].GetCrystal()[1].GetSegID()[0]]));
												EGamma_total = (*Miniball)[0].GetCrystal()[0].GetCore() + (*Miniball)[0].GetCrystal()[1].GetCore();
												EGamDS_total = gamma[idet][ie]*EGamma_total*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[0].GetCluID()][(*Miniball)[0].GetCrystal()[0].GetCryID()][(*Miniball)[0].GetCrystal()[0].GetSegID()[0]]));
												NPT_EGammavsEGamma->Fill((*Miniball)[0].GetCrystal()[0].GetCore()/1.e3,(*Miniball)[0].GetCrystal()[1].GetCore()/1.e3);
												NPT_EDSvsEDS->Fill(EGamDS/1.e3,EGamDS_2/1.e3);
												NPT_EDS_Total->Fill(EGamDS_total/1.e3);
											}
											if ((Miniball->size()==1) && ((*Miniball)[j].GetCrystal().size()==1)){
												NPT_EGamma_Total->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
												NPT_EGamma_1->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
												NPT_EDS_Total->Fill(EGamDS/1.e3);
										 		NPT_EDS_1->Fill(EGamDS/1.e3);
												NPT_EexvsEGamma_1->Fill(Eex[idet][ie]/1.e3,(*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
												NPT_EexvsEDS_1->Fill(Eex[idet][ie]/1.e3,EGamDS/1.e3);
											}
											else if (((Miniball->size()==2) && ((*Miniball)[0].GetCrystal().size()==1) && ((*Miniball)[1].GetCrystal().size()==1)) || ((Miniball->size()==1) && (((*Miniball)[0].GetCrystal().size()==2)))){
												NPT_EGamma_2->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
										 		NPT_EDS_2->Fill(EGamDS/1.e3);
											}
											else{
												NPT_EGamma_3->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
										 		NPT_EDS_3->Fill(EGamDS/1.e3);
											}
											if ((EGamDS<330)&&(310<EGamDS)){
												NPT_EtotvsCmangle_Gamma_320->Fill(theta_cm[idet][ie]/pi*180.,Efor[idet][ie]/1.e3);
												NPT_EtotvsLabangle_Gamma_320->Fill(thetap[idet][ie]/pi*180.,Efor[idet][ie]/1.e3);
												NPT_Eex_Gamma_320->Fill(Eex[idet][ie]/1.e3);
												for (Int_t iAngle=0; iAngle<20.;iAngle++)
												{
													if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_320_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
												}
											}
											if (((EGamDS<340)&&(330<EGamDS))||((EGamDS<310)&&(300<EGamDS))){
												NPT_Eex_Gamma_320_Background->Fill(Eex[idet][ie]/1.e3);
												for (Int_t iAngle=0; iAngle<20.;iAngle++)
												{
													if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_320_Angle_Background[iAngle]->Fill(Eex[idet][ie]/1.e3);
												}
											}
											if (((EGamDS<2160)&&(2060<EGamDS))&&(((*Miniball)[j].GetCrystal()[k].GetCore()<2100)||(2140<(*Miniball)[j].GetCrystal()[k].GetCore()))){
												NPT_EtotvsCmangle_Gamma_2100->Fill(theta_cm[idet][ie]/pi*180.,Efor[idet][ie]/1.e3);
												NPT_Eex_Gamma_2100->Fill(Eex[idet][ie]/1.e3);
												for (Int_t iAngle=0; iAngle<20.;iAngle++)
												{
													if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_2100_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
												}
											}
											if (((EGamDS<2260)&&(2160<EGamDS))&&(((*Miniball)[j].GetCrystal()[k].GetCore()<2100)||(2140<(*Miniball)[j].GetCrystal()[k].GetCore()))){
												NPT_Eex_Gamma_2100_Background->Fill(Eex[idet][ie]/1.e3);
												for (Int_t iAngle=0; iAngle<20.;iAngle++)
												{
													if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_2100_Angle_Background[iAngle]->Fill(Eex[idet][ie]/1.e3);
												}
											}
											if ((EGamDS<2680)&&(2520<EGamDS)){
												NPT_Eex_Gamma_2600->Fill(Eex[idet][ie]/1.e3);
												for (Int_t iAngle=0; iAngle<20.;iAngle++)
												{
													if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_2600_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
												}
											}
											if ((EGamDS<2790)&&(2690<EGamDS)){
												NPT_EtotvsCmangle_Gamma_2700->Fill(theta_cm[idet][ie]/pi*180.,Efor[idet][ie]/1.e3);
												NPT_Eex_Gamma_2700->Fill(Eex[idet][ie]/1.e3);
												for (Int_t iAngle=0; iAngle<20.;iAngle++)
												{
													if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_2700_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
												}
												if (idet==0){
													Q = m[11] + m[Atarget] - m[Apart] - m[11+Atarget-Apart];
													for (Int_t iCD=0;iCD<5;iCD++){
														if (par[2][iCD]==4){
															Pp = sqrt(2*m[Apart]*Efor[idet][ie]);
															Pin = sqrt(2*m[11]*Ein);
															a = 1;
															b = 2*Pp*(cos(thetap[2][j])*cos(thetap[0][ie]) + sin(thetap[2][j])*sin(thetap[0][ie]));
															c = Pp*Pp - Pin*Pin;

															PN = (-b + sqrt(b*b - 4*a*c))/(2*a);
//															theta_coin = (sin(thetap[idet][ie]) - cos(thetap[idet][ie])*tan(thetap[2][iCD]))/(sin(thetain) - cos(thetain)*tan(thetap[2][iCD]));
//															PN = Pp*(theta_coin*cos(thetain) - cos(thetap[idet][ie]))/cos(thetap[2][iCD]);
															EN = pow(PN,2)/(2*m[11+Atarget-Apart]);
//															Pin = Pp*theta_coin;
															Ein_coin = pow(Pin,2)/(2*m[11]);
															Eex_coin = Ein_coin + Q - Efor[idet][ie] - EN;
//														  	Eex_coin = Ein + Q - Efor[idet][ie] - Efor[2][iCD];
															NPT_Eex_coin_Gamma_2700->Fill(Eex_coin/1.e3);
														}
													}	
												}
											}
											if ((EGamDS<3000)&&(2800<EGamDS)){
												NPT_Eex_Gamma_2900->Fill(Eex[idet][ie]/1.e3);
												for (Int_t iAngle=0; iAngle<20.;iAngle++)
												{
													if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_2900_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
												}
											}
											if ((EGamDS<3500)&&(3250<EGamDS)){
												NPT_EtotvsCmangle_Gamma_3300->Fill(theta_cm[idet][ie]/pi*180.,Efor[idet][ie]/1.e3);
												NPT_Eex_Gamma_3300->Fill(Eex[idet][ie]/1.e3);
												for (Int_t iAngle=0; iAngle<20.;iAngle++)
												{
													if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_3300_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
												}
											}
											if ((EGamDS<3500)&&(2500<EGamDS)){
												NPT_Eex_Gamma_3300_comp->Fill(Eex[idet][ie]/1.e3);
												for (Int_t iAngle=0; iAngle<20.;iAngle++)
												{
													if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_3300_comp_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
												}
											}
											if ((EGamDS<3600)&&(3500<EGamDS)){
												NPT_Eex_Gamma_3300_Background->Fill(Eex[idet][ie]/1.e3);
												for (Int_t iAngle=0; iAngle<20.;iAngle++)
												{
													if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_3300_Angle_Background[iAngle]->Fill(Eex[idet][ie]/1.e3);
												}
											}
											if ((EGamDS<6200)&&(4200<EGamDS)){
												NPT_EtotvsCmangle_Gamma_6000->Fill(theta_cm[idet][ie]/pi*180.,Efor[idet][ie]/1.e3);
												NPT_Eex_Gamma_6000_comp->Fill(Eex[idet][ie]/1.e3);
												for (Int_t iAngle=0; iAngle<20.;iAngle++)
												{
													if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_6000_comp_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
												}
											}
											if ((EGamDS<6600)&&(6200<EGamDS)){
												NPT_Eex_Gamma_6000_Background->Fill(Eex[idet][ie]/1.e3);
												for (Int_t iAngle=0; iAngle<20.;iAngle++)
												{
													if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_6000_Angle_Background[iAngle]->Fill(Eex[idet][ie]/1.e3);
												}
											}
										}
									}

									if (strcmp(Histos[ih],"time")==0){
										Gamma_EvsdtNPT->Fill(25*((*Miniball)[j].GetCrystal()[k].GetTime() - time[idet][ie]),(*Miniball)[j].GetCrystal()[k].GetCore()/1000.);
									}
								}

								if (strcmp(Histos[ih],"time")==0){
									// If any particle is seen, fill the total time differences
									if (par[idet][ie]>0){
										Gamma_Evsdttot->Fill(25*((*Miniball)[j].GetCrystal()[k].GetTime() - time[idet][ie]),(*Miniball)[j].GetCrystal()[k].GetCore()/1000.);
										// Cut energies characteristic for inelastic scattering by neutrons on 72Ge, fill time-diff histogram.
										if ((*Miniball)[j].GetCrystal()[k].GetCore()>690.&&(*Miniball)[j].GetCrystal()[k].GetCore()<720.) Gamma_dtGe->Fill(25*((*Miniball)[j].GetCrystal()[k].GetTime() - time[idet][ie]));
										// Cut energies corresponding to another visible line in the total time difference histogram. Beta decay..?
										if ((*Miniball)[j].GetCrystal()[k].GetCore()>2000.&&(*Miniball)[j].GetCrystal()[k].GetCore()<2500.&&(((*Miniball)[j].GetCrystal()[k].GetTime() - time[idet][ie])>0)) Gamma_Ecut->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1000.);
									}
								}
							}
						}
					}
				}
			}
		}
		// End gammas

		if(i%1000 == 0){
			cout<<setw(5)<<setiosflags(ios::fixed)<<setprecision(1)<<(100.*i)/nentries<<" % done\r"<<flush;
		}
	}

	// Writing histograms

	// Single histograms
	for (Int_t ih=0;ih<Histos.size();ih++){
		if (strcmp(Histos[ih],"sipos")==0){
			// Angles
			ForwardBarrel_Theta->Write("",TObject::kOverwrite);
			ForwardBarrel_Phi->Write("",TObject::kOverwrite);
			ForwardBarrel_Theta_det->Write("",TObject::kOverwrite);
			BackwardBarrel_Theta->Write("",TObject::kOverwrite);
			BackwardBarrel_Phi->Write("",TObject::kOverwrite);
			CD_Theta->Write("",TObject::kOverwrite);
			CD_Phi->Write("",TObject::kOverwrite);
			Hitpattern_ForwardBarreltotal->Write("",TObject::kOverwrite);
			Hitpattern_CDtotal->Write("",TObject::kOverwrite);
		}
		
		if (strcmp(Histos[ih],"sie")==0){
			// Energy
			ForwardBarrel_dEvsEtotal->Write("",TObject::kOverwrite);
			BackwardBarrel_dEvsEtotal->Write("",TObject::kOverwrite);
			ForwardBarrel_dEvsEcortotal->Write("",TObject::kOverwrite);
			BackwardBarrel_dEvsEcortotal->Write("",TObject::kOverwrite);
			CD_dEvsEtotal->Write("",TObject::kOverwrite);
			CD_dEvsEcortotal->Write("",TObject::kOverwrite);
		}
		

		if (strcmp(Histos[ih],"part")==0){
			// Particles
			ForwardBarrel_dEvsEcortotal->Write("",TObject::kOverwrite);
			BackwardBarrel_dEvsEcortotal->Write("",TObject::kOverwrite);
			CD_dEvsEcortotal->Write("",TObject::kOverwrite);
			ForwardBarrel_dEvsEcorp->Write("",TObject::kOverwrite);
			ForwardBarrel_dEvsEcord->Write("",TObject::kOverwrite);
			ForwardBarrel_dEvsEcort->Write("",TObject::kOverwrite);
			BackwardBarrel_dEvsEcorp->Write("",TObject::kOverwrite);
			CD_dEvsEcorp->Write("",TObject::kOverwrite);
			CD_dEvsEcord->Write("",TObject::kOverwrite);
			CD_dEvsEcort->Write("",TObject::kOverwrite);
			Proton_EtotvsLabangletotal->Write("",TObject::kOverwrite);
			Proton_Back_EtotvsLabangletotal->Write("",TObject::kOverwrite);
			Proton_EtotvsCmangletotal->Write("",TObject::kOverwrite);
			Deuteron_EtotvsLabangletotal->Write("",TObject::kOverwrite);
			Deuteron_EtotvsCmangletotal->Write("",TObject::kOverwrite);
			Triton_EtotvsLabangletotal->Write("",TObject::kOverwrite);
			Triton_EtotvsCmangletotal->Write("",TObject::kOverwrite);
			Proton_Eextotal->Write("",TObject::kOverwrite);
			Deuteron_Eextotal->Write("",TObject::kOverwrite);
			Triton_Eextotal->Write("",TObject::kOverwrite);
			NPT_Eextotal->Write("",TObject::kOverwrite);
			NPT_Eextotal_Back->Write("",TObject::kOverwrite);
			Proton_EGamma->Write("",TObject::kOverwrite);
			Deuteron_EGamma->Write("",TObject::kOverwrite);
			Triton_EGamma->Write("",TObject::kOverwrite);
			Proton_EGamma_Total->Write("",TObject::kOverwrite);
			Deuteron_EGamma_Total->Write("",TObject::kOverwrite);
			Triton_EGamma_Total->Write("",TObject::kOverwrite);
			Proton_EGamma_1->Write("",TObject::kOverwrite);
			Deuteron_EGamma_1->Write("",TObject::kOverwrite);
			Triton_EGamma_1->Write("",TObject::kOverwrite);
			Proton_EGamma_2->Write("",TObject::kOverwrite);
			Deuteron_EGamma_2->Write("",TObject::kOverwrite);
			Triton_EGamma_2->Write("",TObject::kOverwrite);
			Proton_EGamma_3->Write("",TObject::kOverwrite);
			Deuteron_EGamma_3->Write("",TObject::kOverwrite);
			Triton_EGamma_3->Write("",TObject::kOverwrite);
			Proton_EDS->Write("",TObject::kOverwrite);
			Deuteron_EDS->Write("",TObject::kOverwrite);
			Triton_EDS->Write("",TObject::kOverwrite);
			Proton_EDS_Total->Write("",TObject::kOverwrite);
			Deuteron_EDS_Total->Write("",TObject::kOverwrite);
			Triton_EDS_Total->Write("",TObject::kOverwrite);
			Proton_EDS_1->Write("",TObject::kOverwrite);
			Deuteron_EDS_1->Write("",TObject::kOverwrite);
			Triton_EDS_1->Write("",TObject::kOverwrite);
			Proton_EDS_2->Write("",TObject::kOverwrite);
			Deuteron_EDS_2->Write("",TObject::kOverwrite);
			Triton_EDS_2->Write("",TObject::kOverwrite);
			Proton_EDS_3->Write("",TObject::kOverwrite);
			Deuteron_EDS_3->Write("",TObject::kOverwrite);
			Triton_EDS_3->Write("",TObject::kOverwrite);
			Proton_EexvsEGamma->Write("",TObject::kOverwrite);
			Deuteron_EexvsEGamma->Write("",TObject::kOverwrite);
			Triton_EexvsEGamma->Write("",TObject::kOverwrite);
			Proton_EexvsEGamma_1->Write("",TObject::kOverwrite);
			Deuteron_EexvsEGamma_1->Write("",TObject::kOverwrite);
			Triton_EexvsEGamma_1->Write("",TObject::kOverwrite);
			Proton_EexvsEDS->Write("",TObject::kOverwrite);
			Deuteron_EexvsEDS->Write("",TObject::kOverwrite);
			Triton_EexvsEDS->Write("",TObject::kOverwrite);
			Proton_EexvsEDS_1->Write("",TObject::kOverwrite);
			Deuteron_EexvsEDS_1->Write("",TObject::kOverwrite);
			Triton_EexvsEDS_1->Write("",TObject::kOverwrite);
			Proton_EGammavsEGamma->Write("",TObject::kOverwrite);
			Deuteron_EGammavsEGamma->Write("",TObject::kOverwrite);
			Triton_EGammavsEGamma->Write("",TObject::kOverwrite);
			Proton_EDSvsEDS->Write("",TObject::kOverwrite);
			Deuteron_EDSvsEDS->Write("",TObject::kOverwrite);
			Triton_EDSvsEDS->Write("",TObject::kOverwrite);
			Triton_EDS_coin_3300->Write("",TObject::kOverwrite);
			Deuteron_pr_Time->Write("",TObject::kOverwrite);
		
			// No Punchthrough
			ForwardBarrel_dEvsLabangletotal->Write("",TObject::kOverwrite);
			ForwardBarrel_dEvsCmangletotal->Write("",TObject::kOverwrite);
			BackwardBarrel_dEvsLabangletotal->Write("",TObject::kOverwrite);
			CD_dEvsLabangletotal->Write("",TObject::kOverwrite);
			All_dEvsLabangletotal->Write("",TObject::kOverwrite);
			NPT_EGamma->Write("",TObject::kOverwrite);
			NPT_EGamma_Total->Write("",TObject::kOverwrite);
			NPT_EGamma_1->Write("",TObject::kOverwrite);
			NPT_EGamma_2->Write("",TObject::kOverwrite);
			NPT_EGamma_3->Write("",TObject::kOverwrite);
			NPT_EDS->Write("",TObject::kOverwrite);
			NPT_EDS_Total->Write("",TObject::kOverwrite);
			NPT_EDS_1->Write("",TObject::kOverwrite);
			NPT_EDS_2->Write("",TObject::kOverwrite);
			NPT_EDS_3->Write("",TObject::kOverwrite);
			NPT_EGamma_Low->Write("",TObject::kOverwrite);
			NPT_EDS_Low->Write("",TObject::kOverwrite);
			NPT_EGamma_High->Write("",TObject::kOverwrite);
			NPT_EDS_High->Write("",TObject::kOverwrite);
			NPT_EGamma_Back->Write("",TObject::kOverwrite);
			NPT_EDS_Back->Write("",TObject::kOverwrite);
			NPT_EexvsEGamma->Write("",TObject::kOverwrite);
			NPT_EexvsEDS->Write("",TObject::kOverwrite);
			NPT_EexvsEGamma_Back->Write("",TObject::kOverwrite);
			NPT_EexvsEDS_Back->Write("",TObject::kOverwrite);
			NPT_EexvsEGamma_1->Write("",TObject::kOverwrite);
			NPT_EexvsEDS_1->Write("",TObject::kOverwrite);
			NPT_EGammavsEGamma->Write("",TObject::kOverwrite);
			NPT_EDSvsEDS->Write("",TObject::kOverwrite);

			// 11Be
			Be11_Angle->Write("",TObject::kOverwrite);
			Be11_pr_Time->Write("",TObject::kOverwrite);
			
			Gamma_EDS->Write("",TObject::kOverwrite);
		
			//Flour_Eextotal->Write("",TObject::kOverwrite);
			
			Proton_Eex_Gamma_511->Write("",TObject::kOverwrite);
			Proton_Eex_Gamma_511_Background->Write("",TObject::kOverwrite);
			Proton_Eex_Gamma_2100->Write("",TObject::kOverwrite);
			Proton_Eex_Gamma_2100_1->Write("",TObject::kOverwrite);
			Proton_Eex_Gamma_2700->Write("",TObject::kOverwrite);
			Proton_Eex_Gamma_2700_1->Write("",TObject::kOverwrite);
			Deuteron_Eex_Gamma_320->Write("",TObject::kOverwrite);
			Deuteron_Eex_Gamma_320_1->Write("",TObject::kOverwrite);
			Triton_Eex_Gamma_3300->Write("",TObject::kOverwrite);
			Triton_Eex_Gamma_3300_comp->Write("",TObject::kOverwrite);
			Triton_Eex_Gamma_3300_1->Write("",TObject::kOverwrite);
			Triton_Eex_Gamma_6000_comp->Write("",TObject::kOverwrite);
			Triton_Eex_Gamma_6000_comp_1->Write("",TObject::kOverwrite);
			NPT_Eex_Gamma_320->Write("",TObject::kOverwrite);
			NPT_Eex_Gamma_320_Background->Write("",TObject::kOverwrite);
			NPT_Eex_Gamma_2100->Write("",TObject::kOverwrite);
			NPT_Eex_Gamma_2100_Background->Write("",TObject::kOverwrite);
			NPT_Eex_Gamma_2600->Write("",TObject::kOverwrite);
			NPT_Eex_Gamma_2700->Write("",TObject::kOverwrite);
			NPT_Eex_Gamma_2900->Write("",TObject::kOverwrite);
			NPT_Eex_Gamma_3300->Write("",TObject::kOverwrite);
			NPT_Eex_Gamma_3300_comp->Write("",TObject::kOverwrite);
			NPT_Eex_Gamma_3300_Background->Write("",TObject::kOverwrite);
			NPT_Eex_Gamma_6000_comp->Write("",TObject::kOverwrite);
			NPT_Eex_Gamma_6000_Background->Write("",TObject::kOverwrite);
			
			NPT_Eex_coin_Gamma_2700->Write("",TObject::kOverwrite);
				
			Proton_EtotvsCmangle_Gamma_2100->Write("",TObject::kOverwrite);
			Proton_EtotvsCmangle_Gamma_2700->Write("",TObject::kOverwrite);
			Deuteron_EtotvsCmangle_Gamma_320->Write("",TObject::kOverwrite);
			Triton_EtotvsCmangle_Gamma_3300->Write("",TObject::kOverwrite);
			Triton_EtotvsCmangle_Gamma_6000->Write("",TObject::kOverwrite);
			NPT_EtotvsCmangle_Gamma_2100->Write("",TObject::kOverwrite);
			NPT_EtotvsCmangle_Gamma_2700->Write("",TObject::kOverwrite);
			NPT_EtotvsCmangle_Gamma_320->Write("",TObject::kOverwrite);
			NPT_EtotvsLabangle_Gamma_320->Write("",TObject::kOverwrite);
			NPT_EtotvsCmangle_Gamma_3300->Write("",TObject::kOverwrite);
			NPT_EtotvsCmangle_Gamma_6000->Write("",TObject::kOverwrite);
			
			Proton_MBTheta_2100->Write("",TObject::kOverwrite);
			Proton_MBTheta_2700->Write("",TObject::kOverwrite);
			Deuteron_MBTheta_320->Write("",TObject::kOverwrite);
			Triton_MBTheta_3300->Write("",TObject::kOverwrite);
		
	//		Proton_EtotvsLabangle_Gamma_1000->Write("",TObject::kOverwrite);
		}

		if (strcmp(Histos[ih],"coin")==0){
			Hitpattern_cointotal->Write("",TObject::kOverwrite);
			XY_cointotal->Write("",TObject::kOverwrite);
			XY_Beam_cointotal->Write("",TObject::kOverwrite);
			DetvsDet_coin->Write("",TObject::kOverwrite);
			Proton_Eex_cointotal->Write("",TObject::kOverwrite);
			Deuteron_Eex_cointotal->Write("",TObject::kOverwrite);
			Triton_Eex_cointotal->Write("",TObject::kOverwrite);
			NPT_Eex_cointotal->Write("",TObject::kOverwrite);
			Triton_Ein1_cointotal->Write("",TObject::kOverwrite);
			Triton_Ein2_cointotal->Write("",TObject::kOverwrite);
			NPT_EexvsEex1_cointotal->Write("",TObject::kOverwrite);
			NPT_EexvsEex2_cointotal->Write("",TObject::kOverwrite);
			NPT_EtotvsLabangletotal_coin->Write("",TObject::kOverwrite);


		}	
	

	
	
		if (strcmp(Histos[ih],"gee")==0){
			// Gamma histograms
			Gamma_Etotal->Write("",TObject::kOverwrite);
		}
	
		// Pad dependent histograms
		for (Int_t i=0;i<4;i++){
			
			if (strcmp(Histos[ih],"sipos")==0){
				// Position
				ForwardBarrel_PosvsStrip[i]->Write("",TObject::kOverwrite);
				BackwardBarrel_PosvsStrip[i]->Write("",TObject::kOverwrite);
			}

			if (strcmp(Histos[ih],"sie")==0){
				// dE
				ForwardBarrel_dEtot[i]->Write("",TObject::kOverwrite);
				BackwardBarrel_dEtot[i]->Write("",TObject::kOverwrite);
				ForwardBarrel_dEvsStrip[i]->Write("",TObject::kOverwrite);
				BackwardBarrel_dEvsStrip[i]->Write("",TObject::kOverwrite);
				BackwardBarrel_dEvsPos[i]->Write("",TObject::kOverwrite);
				ForwardBarrel_dEvsDetangle[i]->Write("",TObject::kOverwrite);
		
				// E backdetector
				ForwardBarrel_Eback[i]->Write("",TObject::kOverwrite);
				BackwardBarrel_Eback[i]->Write("",TObject::kOverwrite);
				ForwardBarrel_dEvsEtot[i]->Write("",TObject::kOverwrite);
				BackwardBarrel_dEvsEtot[i]->Write("",TObject::kOverwrite);
				ForwardBarrel_dEvsEcortot[i]->Write("",TObject::kOverwrite);
				BackwardBarrel_dEvsEcortot[i]->Write("",TObject::kOverwrite);
			}

			if (strcmp(Histos[ih],"alpha")==0){
				// Efor
				ForwardBarrel_Efor[i]->Write("",TObject::kOverwrite);
				ForwardBarrel_EforvsDetangle[i]->Write("",TObject::kOverwrite);
			}

			if (strcmp(Histos[ih],"part")==0){
				ForwardBarrel_dEvsEcortot[i]->Write("",TObject::kOverwrite);
				BackwardBarrel_dEvsEcortot[i]->Write("",TObject::kOverwrite);
				
				// No punchthrough
				ForwardBarrel_dEvsLabangle[i]->Write("",TObject::kOverwrite);
				BackwardBarrel_dEvsLabangle[i]->Write("",TObject::kOverwrite);
		
				// Particles
				Proton_EtotvsLabangle[i]->Write("",TObject::kOverwrite);
				Deuteron_EtotvsLabangle[i]->Write("",TObject::kOverwrite);
				Triton_EtotvsLabangle[i]->Write("",TObject::kOverwrite);
				Proton_Eex_ForwardBarrel[i]->Write("",TObject::kOverwrite);
				Deuteron_Eex_ForwardBarrel[i]->Write("",TObject::kOverwrite);
				Triton_Eex_ForwardBarrel[i]->Write("",TObject::kOverwrite);
				Proton_Eex_CD[i]->Write("",TObject::kOverwrite);
				Deuteron_Eex_CD[i]->Write("",TObject::kOverwrite);
				Triton_Eex_CD[i]->Write("",TObject::kOverwrite);
				NPT_Eex_BackwardBarrel[i]->Write("",TObject::kOverwrite);
				Proton_EexvsLabangle_ForwardBarrel[i]->Write("",TObject::kOverwrite);
				Deuteron_EexvsLabangle_ForwardBarrel[i]->Write("",TObject::kOverwrite);
				Triton_EexvsLabangle_ForwardBarrel[i]->Write("",TObject::kOverwrite);
				Proton_EexvsLabangle_CD[i]->Write("",TObject::kOverwrite);
				Deuteron_EexvsLabangle_CD[i]->Write("",TObject::kOverwrite);
				Triton_EexvsLabangle_CD[i]->Write("",TObject::kOverwrite);
			}
			
			if (strcmp(Histos[ih],"coin")==0){
				// Coincidence
				Proton_Eex_coin[i]->Write("",TObject::kOverwrite);
				Deuteron_Eex_coin[i]->Write("",TObject::kOverwrite);
				Triton_Eex_coin[i]->Write("",TObject::kOverwrite);
				Proton_EexvsLabangle_coin[i]->Write("",TObject::kOverwrite);
				Deuteron_EexvsLabangle_coin[i]->Write("",TObject::kOverwrite);
				Triton_EexvsLabangle_coin[i]->Write("",TObject::kOverwrite);
				Proton_ThetaIn_coin[i]->Write("",TObject::kOverwrite);
				Deuteron_ThetaIn_coin[i]->Write("",TObject::kOverwrite);
				Triton_ThetaIn_coin[i]->Write("",TObject::kOverwrite);
				Triton_PhiIn_coin[i]->Write("",TObject::kOverwrite);
				Proton_ThetaInvsEex_coin[i]->Write("",TObject::kOverwrite);
				Deuteron_ThetaInvsEex_coin[i]->Write("",TObject::kOverwrite);
				Triton_ThetaInvsEex_coin[i]->Write("",TObject::kOverwrite);
				Proton_PhivsPhi_coin[i]->Write("",TObject::kOverwrite);
				Deuteron_PhivsPhi_coin[i]->Write("",TObject::kOverwrite);
				Triton_PhivsPhi_coin[i]->Write("",TObject::kOverwrite);
		
				Hitpattern_coin[i]->Write("",TObject::kOverwrite);
				XY_coin[i]->Write("",TObject::kOverwrite);
				XY_Beam_coin[i]->Write("",TObject::kOverwrite);
			}	
	
			// Strip dependent histograms
			for (Int_t j=0;j<16;j++){
				if (strcmp(Histos[ih],"sipos")==0){
					// Position
					ForwardBarrel_Pos[i][j]->Write("",TObject::kOverwrite);
					BackwardBarrel_Pos[i][j]->Write("",TObject::kOverwrite);
				}

				if (strcmp(Histos[ih],"sie")==0){
					// dE
					ForwardBarrel_dE[i][j]->Write("",TObject::kOverwrite);
					ForwardBarrel_dEvsPos[i][j]->Write("",TObject::kOverwrite);
					for (Int_t k=0;k<16;k++) ForwardBarrel_dEforPos[i][j][k]->Write("",TObject::kOverwrite);
					BackwardBarrel_dE[i][j]->Write("",TObject::kOverwrite);
					BackwardBarrel_dEforPos[i][j]->Write("",TObject::kOverwrite);
		
					// E backdetector
					ForwardBarrel_dEvsE[i][j]->Write("",TObject::kOverwrite);
					ForwardBarrel_dEvsEcor[i][j]->Write("",TObject::kOverwrite);
					BackwardBarrel_dEvsE[i][j]->Write("",TObject::kOverwrite);
				}


				if (strcmp(Histos[ih],"alpha")==0){
					// Efor
					ForwardBarrel_EforvsPos[i][j]->Write("",TObject::kOverwrite);
				}
			}
	
			if (strcmp(Histos[ih],"sie")==0){
				// Angle dependent histograms
				for (Int_t j=0;j<20;j++){
					// dE
					ForwardBarrel_dEangle[i][j]->Write("",TObject::kOverwrite);
					ForwardBarrel_dEvsEangle[i][j]->Write("",TObject::kOverwrite);
					BackwardBarrel_dEvsEangle[i][j]->Write("",TObject::kOverwrite);
				}
			}
		}
	
		if (strcmp(Histos[ih],"part")==0){
			// Angle dependent histograms
			for (Int_t i=0;i<20;i++){
					Proton_Eex_Angle[i]->Write("",TObject::kOverwrite);
					Proton_Eex_511_Angle[i]->Write("",TObject::kOverwrite);
					Proton_Eex_511_Angle_Background[i]->Write("",TObject::kOverwrite);
					Proton_Eex_2100_Angle[i]->Write("",TObject::kOverwrite);
					Proton_Eex_2100_Angle_Background[i]->Write("",TObject::kOverwrite);
					Proton_Eex_2700_Angle[i]->Write("",TObject::kOverwrite);
					Deuteron_Eex_Angle[i]->Write("",TObject::kOverwrite);
					Deuteron_Eex_320_Angle[i]->Write("",TObject::kOverwrite);
					Deuteron_Eex_320_Angle_Background[i]->Write("",TObject::kOverwrite);
					Triton_Eex_Angle[i]->Write("",TObject::kOverwrite);
					Triton_Eex_3300_Angle[i]->Write("",TObject::kOverwrite);
					Triton_Eex_3300_comp_Angle[i]->Write("",TObject::kOverwrite);
					Triton_Eex_2600_Angle[i]->Write("",TObject::kOverwrite);
					Triton_Eex_2900_Angle[i]->Write("",TObject::kOverwrite);
					Triton_Eex_6000_Angle[i]->Write("",TObject::kOverwrite);
					Triton_Eex_6000_comp_Angle[i]->Write("",TObject::kOverwrite);
					NPT_Eex_Angle[i]->Write("",TObject::kOverwrite);
					NPT_Eex_Angle_Back[i]->Write("",TObject::kOverwrite);
					NPT_Eex_320_Angle[i]->Write("",TObject::kOverwrite);
					NPT_Eex_320_Angle_Background[i]->Write("",TObject::kOverwrite);
					NPT_Eex_2100_Angle[i]->Write("",TObject::kOverwrite);
					NPT_Eex_2100_Angle_Background[i]->Write("",TObject::kOverwrite);
					NPT_Eex_2600_Angle[i]->Write("",TObject::kOverwrite);
					NPT_Eex_2700_Angle[i]->Write("",TObject::kOverwrite);
					NPT_Eex_2900_Angle[i]->Write("",TObject::kOverwrite);
					NPT_Eex_3300_Angle[i]->Write("",TObject::kOverwrite);
					NPT_Eex_3300_comp_Angle[i]->Write("",TObject::kOverwrite);
					NPT_Eex_3300_Angle_Background[i]->Write("",TObject::kOverwrite);
					NPT_Eex_6000_comp_Angle[i]->Write("",TObject::kOverwrite);
					NPT_Eex_6000_Angle_Background[i]->Write("",TObject::kOverwrite);
					Beryllium_Eex_Angle_Stop_FB[i]->Write("",TObject::kOverwrite);
					Beryllium_Eex_Angle_Stop_CD[i]->Write("",TObject::kOverwrite);
			}
		}
		
	
		// Cd dependent histograms	
		for (Int_t i=0;i<4;i++){
				
			if (strcmp(Histos[ih],"sipos")==0){
				// pos
				CD_StripvsRing[i]->Write("",TObject::kOverwrite);
			}
			if (strcmp(Histos[ih],"sie")==0){
				// dE
				CD_dEtot[i]->Write("",TObject::kOverwrite);
				CD_dEvsStrip[i]->Write("",TObject::kOverwrite);
				CD_dEvsRing[i]->Write("",TObject::kOverwrite);
				CD_EStripvsERing[i]->Write("",TObject::kOverwrite);
				CD_dEvsEtot[i]->Write("",TObject::kOverwrite);
				CD_dEvsEcortot[i]->Write("",TObject::kOverwrite);
				CD_Eback[i]->Write("",TObject::kOverwrite);
			}

			if (strcmp(Histos[ih],"part")==0){
				// Angle
				CD_dEvsLabangle[i]->Write("",TObject::kOverwrite);
	
			}
	
			if (strcmp(Histos[ih],"sie")==0){
				// Strip/ring dependent histograms
				for (Int_t j=0;j<16;j++){
					CD_dEforStrip[i][j]->Write("",TObject::kOverwrite);
					CD_dEforRing[i][j]->Write("",TObject::kOverwrite);
					CD_dEvsEStrip[i][j]->Write("",TObject::kOverwrite);
					CD_dEvsERing[i][j]->Write("",TObject::kOverwrite);
				}
			}
		}
	
	
	
		if (strcmp(Histos[ih],"time")==0){
			// Time
			Gamma_EvsdtProton->Write("",TObject::kOverwrite);
			Proton_dt_511->Write("",TObject::kOverwrite);
			Gamma_EvsdtDeuteron->Write("",TObject::kOverwrite);
			Gamma_EvsdtTriton->Write("",TObject::kOverwrite);
			Gamma_EvsdtNPT->Write("",TObject::kOverwrite);
			Gamma_EvsdtBack->Write("",TObject::kOverwrite);
			Gamma_EvsdtCD->Write("",TObject::kOverwrite);
			Gamma_Evsdttot->Write("",TObject::kOverwrite);
			Gamma_dtGe->Write("",TObject::kOverwrite);
			Gamma_Ecut->Write("",TObject::kOverwrite);
		}

		// Clusters
		for (Int_t i=0;i<8;i++){
		
			// Crystals
			for (Int_t j=0;j<3;j++){
				if (strcmp(Histos[ih],"gee")==0){
					Gamma_ECrystals1[i][j]->Write("",TObject::kOverwrite);
					Gamma_ECrystals2[i][j]->Write("",TObject::kOverwrite);
				}

				// Segments
				for (Int_t k=0;k<6;k++){
					if (strcmp(Histos[ih],"gee")==0) Gamma_ESegments[i][j][k]->Write("",TObject::kOverwrite);
					if (strcmp(Histos[ih],"part")==0){
						Proton_ESegments[i][j][k]->Write("",TObject::kOverwrite);
						Deuteron_ESegments[i][j][k]->Write("",TObject::kOverwrite);
					}
				}
			}
		}
	}
	outfile->Close();
	delete tr;


	printf("Total: Np = %d, Nd = %d, Nt = %d, t = %f min \n",Np/Histos.size(),Nd/Histos.size(),Nt/Histos.size(),EbisTime*25.e-9/60.);
	printf("Pr. Minute: Np = %f, Nd = %f, Nt = %f \n",Np/Histos.size()/(EbisTime*25.e-9/60.),Nd/Histos.size()/(EbisTime*25.e-9/60.),Nt/Histos.size()/(EbisTime*25.e-9/60.));
	printf("Pr. Minute in FB: Np = %f, Nd = %f, Nt = %f \n",Np_FB/Histos.size()/(EbisTime*25.e-9/60.),Nd_FB/Histos.size()/(EbisTime*25.e-9/60.),Nt_FB/Histos.size()/(EbisTime*25.e-9/60.));
	printf("Pr. Minute in CD: Np = %f, Nd = %f, Nt = %f \n",Np_CD/Histos.size()/(EbisTime*25.e-9/60.),Nd_CD/Histos.size()/(EbisTime*25.e-9/60.),Nt_CD/Histos.size()/(EbisTime*25.e-9/60.));


	Nfile.open("Npar.txt", fstream::out | fstream::app);

	Nfile << Np/(EbisTime*25.e-9/60.) << "	" << Nd/(EbisTime*25.e-9/60.) << "	" << Nt/(EbisTime*25.e-9/60.) << endl;

	Nfile.close();
	
	Nfile.open("Npar_FB.txt", fstream::out | fstream::app);

	Nfile << Np_FB/(EbisTime*25.e-9) << "	" << Nd_FB/(EbisTime*25.e-9) << "	" << Nt_FB/(EbisTime*25.e-9) << endl;

	Nfile.close();
	
	Nfile.open("Npar_CD.txt", fstream::out | fstream::app);

	Nfile << Np_CD/(EbisTime*25.e-9/60.) << "	" << Nd_CD/(EbisTime*25.e-9/60.) << "	" << Nt_CD/(EbisTime*25.e-9/60.) << endl;

	Nfile.close();
	
	Nfile.open("Time.txt", fstream::out | fstream::app);

	Nfile << EbisTime*25.e-9/60. << endl;

	Nfile.close();
	
	return 0;

}

